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ABSTRACT 


A review and summary of recent developments of the Green's function 
method and the computer program SOUSSA (Steady, Oscillatory, and ynsteady Subsonic 
and Supersonic Aerodynamics) are presented. Applying the Green's function method to 
the fully unsteady (transient) potential equation yields an integro-differentia! -delay 
equation. With spatial discretization by the finite-element method, this equation is 
approximated by a set of differential -delay equations in time. Time solution by 
Laplace transform yields a matrix relating the velocity potential to the normal wash. 
Premultiplying and postmultiplying by the matrices relating generalized forces to the 
potential and the normal wash to the generalized coordinates one obtains the matrix 
of the generalized aerodynamic forces. The frequency and mode-shape dependence of 
this matrix makes the program SOUSSA very useful for multiple frequency and repeated 
mode-shape evaluations. The program SOUSSA is general, flexible, easy to use, and 
accurate. Applications to aerodynamic design are also discussed. The user/program- 
mer manual for SOUSSA-P 1 . 1 is presented in Volume 2 of this report (NASA CR- 
159131). 
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SECTION 1 


INTRODUCTION 

Presented here is a general formulation for steady, oscillatory, and unsteady 
subsonic and supersonic, potential aerodynamics for complex aircraft configurations. 
This formulation is the basis for the computer program SOUSSA-P, which may be used 
for a variety of aerodynamic computations. These include, for example, the follow- 
ing types of aerodynamic analysis: 

• Unsteady State Applications 

a) Flutter or Gust Analysis. 

b) Flutter or Gust Analysis with multiple sets of frequencies. 

c) Flutter or Gust Analysis with multiple sets of boundary- 
condition modes and/or generalized -force modes. 

d) Flutter or Gust Analysis with multiple Mach numbers. 

• Steady or Quasi -Steady State Applications 

a) Steady-State Pressure Distributions. 

b) Structural Design Loads. 

c) Aerodynamic Coefficients. 

d) Stability Derivatives. 

e) Static Aeroelastic Analysis. 

Reviews of other methods in this field are given in Refs. 1 and 2 and, there- 
fore, are not presented here. The method is based upon the Green's function method 
(to transform the velocity-potential differential equation into an integral-differential- 
delay equation) and the finite element method (to transform the equation into a set of 
differential -delay equations in time). 
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A review of the development of the method is presented in the following 
subsection. The integral equation used in SOUSSA-P is obtained in Section 2, whereas 
the numerical formulation for the solution of the integral equation is presented in 
Section 3. The boundary conditions are considered in Section 4, whereas Section 5 
deals with the evaluation of pressure coefficients from the potential and the evaluation 
of the generalized forces from the pressure coefficient. Section 6 discusses the matrix 
relating the generalized aerodynamic forces to the generalized coordinates describing 
the rigid-body motion and/or deformation of the aircraft. The use of the method in 
aerodynamic design (find the shape given the pressure distribution) is also discussed in 
Section 6. Concluding remarks are presented in Section 7. 

Subtle points in the formulation are dealt with in the Appendices. In parti- 
cular, the value of the function E for a point on the surface is obtained in Appendix 
A. A discussion of Laplace transform and the truncation of the wake is given in Appen- 
dix B. The effect of the motion of the surface of the body is discussed in Appendix 
C. Appendix D deals with the boundary conditions for a body-axis formulation. In 
Appendix E the closed-wake phenomenon is examined. 

1 . 1 REVIEW OF DEVELOPMENT OF METHOD 

The method is based upon a formulation developed by Morino (Refs. 3 and 4) 
in which the Green’s function method is applied to the equation for the velocity poten- 
tial. Use of the infinite-space Green's function method yields a representation of the 
potential cp at any point, p , in the flow field (control point) in terms of the values of 
the potential and its normal derivatives on the surface, E, surrounding the body and its 
wake. The integral equation is obtained by imposing the requirement that the value of 
the potential at p approaches the value of cp on the surface if p approaches a point on 
the surface. The wake is a natural byproduct of the method and is treated as a layer of 
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doublets. It may be noted that the Integral equation is different from other ones cur- 
rently used in aerodynamics in that it does not require that the boundary condition on 
the normal wash be satisfied, but rather makes use of the continuity of the potential as 
the control point approaches the surface, a. The two boundary conditions are mathe- 
matical equivalents, however. The tangency boundary conditions are automatically 
satisfied by the type of representation obtained with the Green's theorem. A feasibility 
analysis by Morino and Kuo (Refs. 5 and 6) indicates that the method is fast and accu- 
rate. A general formulation for steady and oscillatory subsonic and supersonic flow is 
given in Ref. 7. In that formulation the surface of the aircraft and its wake was 
divided into small quadrilateral elements. Each element was replaced by a portion of 
a hyperboloidal paraboloid (hyperboloidal element) defined by the four corner points of 
the element. In this process, the continuity of the surface was maintained, although 
discontinuities in the slopes were introduced. The unknown potential was assumed to 
be constant within each element (zeroth-order element) and therefore, the integral 
equation was approximated with a system of algebraic equations. 

An extension to the fully unsteady formulation, both in the time and com- 
plex-frequency (Laplace) domain, is presented in Ref. 8 for an arbitrary finite-element 
representation. Evaluation of the matrix of the generalized aerodynamic forces is pre- 
sented in Ref. 9. Modification of the theory for coplanar surfaces is discussed in Ref. 
10. Numerical results for subsonic and supersonic flows have been obtained for wing- 
body-tail configurations in steady oscillatory and full unsteady (transient) subsonic and 
supersonic flows (Refs. 9 and 10). 

More recently the computer time for the evaluation of the aerodynamic in- 
fluence coefficients has been reduced considerably by using numerical quadrature to 
evaluate the source and doublet integrals for "distant" elements (Ref. 11). Additional 
work includes the wake roll-up (Ref. 12). An overall review of the method is given in 
Ref. 13. An assessment of the method is given in Ref. 14. 
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1.2 THE COMPUTER PROGRAM SOUSSA-P 1 . 1 

The method reviewed in Section 1.1 was originally imbedded in the general- 
purpose computer program SOUSSA-I (Steady, Oscillatory, and Unsteady Subsonic and 
Supersonic Aerodynamics - Interim Version), which was first presented in Referece 9. 
SOUSSA-1 was at the time the most general code available for potential aerodynamic 
analysis. However, SOUSSA-I was developed in an academic environment, for research 
and proof-of-concept purposes. 'A more general and efficient version (SOUSSA-P) was 
planned for "production" applications. The SOUSSA-P program possesses the following 
qualitites: 

• User Orientation - To enable its use without extensive specialized 
training. 

• Generality - The program has been structured to facilitate the 
analysis of aerodynamic problems involving a wide range of flight 
speeds, arbitrary aircraft geometry, multiple sets of vibration and 
deformation modes, and multiple sets of frequencies. Furthermore, 
SOUSSA-P was designed to be compatible with most currently avail- 
able geometry preprocessors. 

• Computational Efficiency - A premium has been placed on the con- 
servation of central processor time as well as central memory so that 
the program may be a useful tool for application to complicated con- 
figurations such as complete aircraft. 

• Simplicity of Method - It provides a commonality between subsonic 
and supersonic, steady and unsteady flows. Also, the expressions for 
the coefficients are very simple due to the vector formulation of the 
problem. Furthermore, the use of quadrilateral hyperboloidal elements, 
described in terms of their corner points, allows for application to 
arbitrarily complex configurations. 

• Accuracy - The method is accurate and fast, as indicated by the 
results obtained thus far. 

• Modularity - To facilitate incorporating new or improved features of 
the method or the computational algorithms. 

SOUSSA-P employs the sophisticated data handling capabilities of the SPAR 
Finite Element Structural Analysis, System Level II, computer program (Ref. 15), which 
is in part responsible for realizing the aforementioned qualities. 
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1.3 


FORMULATION OF PROBLEM 


The integral equation used in SOUSSA-P is derived in Ref. 3. For the 
sake of clarity, a simpler derivation is presented in this report, assuming from the 
beginning that the motion of the surface is infinitesimal.* 

In this report, the isentropic inviscid flow of a perfect gas, initially 
irrotational , is considered. Under this hypothesis, the flow can be described by the 
velocity potential, 0 , such that the velocity Vp of the fluid is given by 



0 - 1 ) 


where 7 is the del operator (nabla). Consider a frame of reference i, j, k such that the 
undisturbed flow has velocity U in the direction of the i axis. Then the small - 

' CD 

perturbation, linearized equation of the unsteady aerodynamic potential is 




0 - 2 ) 


2 

where Cj^ is the speed of sound of the undisturbed flow, v is the Laplacian operator, 
and 


_d_ 

dt 


b_ 

at 


+ u 

00 


a 

a x 


0-3) 


is the linear substantial derivative. Then, it is convenient to introduce the perturbation 
potential 9 , such that 

0 = U^x +9 (1-4) 

Note that 

cp=0 (1-5) 

in the undisturbed flow. 

*This point is discussed further in Appendix C. 
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A very general approach is considered in Ref. 3 by assuming that the body 
immersed in this flow has arbitrary shape and is moving with arbitrary motion, including 
deformation. Thus, the surface of the body is represented in the general form 

S (x, y, z, t) = 0 (1 -6) 


The no -flow -through boundary condition on the body is given by (see 
Subsection 4-1) 


BS 

— + 70 . VS = 0 

a t 


on S = 0 


(1-7) 


By using Equations (1-3) and (1-4), Equation (1-7) reduces to* 


7 <p.VS = -f — +IL 


at 


oo 


ax 


(1-8) 


or 


*B*e. = — 1 

3n | 9 s | 


dS 

dt 


(1-9) 


As mentioned previously, a simplified approach is used here, in which the surface is 
assumed to move only infinitesimally with respect to the frame of reference traveling at 
velocity (in direction of the negative x-axis) with respect to the undisturbed flow. 
Mathematically, it is assumed that the surface of the body is represented in the form 

S (x, y, z, t) = S B (x, y, z) + eS^x, y, z, t) =0 (1-10) 


*This is the exact boundary condition: no small -perturbation assumption is used here, 
A nonzero flow-through boundary condition can be easily included in the formulation 
by adding the flow-through term in Equation (1-8). 
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with e < < 1 . In Equation (1-10), Sg represents the steady-state geometry while Sy 
gives the unsteady perturbation. 

Next consider the pressure distribution which can be evaluated from the 
exact nonlinear Bernoulli Theorem for barotropic fluids 

f JL dp = -I OH + U ?£ + 1 Vcp • Vcp) (1-11) 

Poo P \3t 00 3x 2 j 

or if nonlinear terms are neglected. 


C . JL2=. 
P 1 p U 2 

2 oo 


= -2 


ir 


^ + u_ 

i3t 


oo 


3cp 

bxi 


2 dcp 


oo 


u 


2 dt 
'oo 


^ (1-12) 


which is the form used in SOUSSA-P 1.1. 

Then the nondimensional generalized aerodynamic force acting on the non- 
dimensional mode Mg is given by 


■4 0 


C n • M. da R 
o J? P h B 

* CT B 


0-13) 


where <jg is the surface of the body and l is a reference length. 

Finally, consider the boundary condition on the wake, i.e., that the pressure 
discontinuity across the wake be equal to zero. Using Equation (1-11) one obtains* 



( v ?y + Vcp^) • V 


(?u - Vj) 


0 


(1-14) 


*The subscripts u and i, indicate the upper and lower sides of the wake. They are 
used even when the wake is vertical or rolled-up. 
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or 


Acp = 0 


(1-15) 


_D 

d7 

where Acp = cp o - cp^ and D/Dt is the substantial derivative, i.e., the time deriva- 
tive evaluated by following a wake particle which by definition has velocity 

V = T + 1 rv Tu + v % ) (H6) 

Note that for steady state Equation (1-15) becomes 

_1_ Acp = 0 (1-17) 

ds W 

where s^y is the arc length along a wake streamline. Equations (1-15) and (1-17) 
are exact, i.e., involve no small-perturbation assumption. Equation (1-15) may be 
integrated to yield 

Acp(p, t) = Acp(p TE , t - t) ( 1 - 18 ) 

where t is the time necessary for a particle in the wake to travel from the point, p-j-^ 
(origin of the vortex line at the trailing edge), to the point p . Small perturbation 
hypothesis is used in SOUSSA-P 1.1 and hence t is approximated by 

T = (X - X^)/^ (1-19) 

This is consistent with using the linearized Bernoulli's theorem for the evaluation of 
the pressure since in this case Equation (1-12) yields 
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or 


Acp - A cp(x - t) 


( 1 - 20 ) 


( 1 - 21 ) 


An additional boundary condition is that the normal component of the 
velocity is zero on the wake, i.e., 

= 0 ( 1 - 22 ) 

dn 

This condition may be used to obtain the geometry of the wake. 
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SECTION 2 


SUBSONIC AND SUPERSONIC INTEGRAL EQUATIONS 

Considered in this section is the integral equation for unsteady subsonic 
and supersonic potential aerodynamics for an aircraft having arbitrary shape. The rigid- 
body motion and/or deformation of the aircraft is assumed to consist of small perturba- 
tions (starting at t = 0) with respect to the constant-speed motion. The objective of 
this formulation is to describe the functional relationship between aerodynamic poten- 
tial and its normal derivative (normal wash, \|f = dcp/dn) on the fluid boundary. 

The analysis presented in this section is based upon an integral formulation 
presented in Refs. 3 and 4, which includes completely arbitrary motion. Note that 
in order to perform a linear-system analysis of the aircraft, it is convenient to use a 
general aerodynamic formulation, i.e., fully transient response for time-domain 
analysis, and the aerodynamic transfer function (Laplace transform of the fully unsteady 
operator) for complex-frequency-domain analysis. Consistent with this type of analysis, 
the unsteady contribution is assumed to start at time t = 0, so that for time t < 0 the 
flow is in steady state. Furthermore, the unsteadiness of the aircraft is assumed to 
consist of small (infinitesimal) perturbations around the steady-state configuration. 

In this section the subsonic formulation is presented in detail; the extension 
to supersonic flows is also outlined. 

2. 1 GREEN 'S THEOREM FOR POTENTIAL AERODYNAMICS 

The purpose of this analysis is to obtain a representation of the potential in 
terms of its value (and the values of its derivatives) on the surface of the body and the 
wake using Green's Theorem. 
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Note that the equation of the aerodynamic potential given by Equation (1-2) 
is not valid on the wake, where discontinuities on cp exist. Thus, consider the volume V 
in which Equation (1-2) is valid. At any instant of time, this volume is given by the 
whole physical space except the volume, Vg, occupied by the body and the infinites- 
imally thin layer, Vyy, representing the wake. As mentioned above, the volume V is 
assumed to be time-independent. Define the function E (see Figure 2-1). 

E(x, y, z) = 1 in V 

= 0 mV B + V w (2-1) 

This function represents the domain of validity of the equation of the potential and will 
be called "domain function." Consider the surface of discontinuity of the function E, 
that is, the surface, cr, surrounding the volume Vg + . Let 

tyx, y, z)=° (2-2) 

be the equation of the surface o. 

Note that the surface ct is composed of two branches. The first, ag, is the 
surface of the body given by Equation (1 -10) (with e = 0 except for the boundary con- 
ditions*), The second is the surface, Oy^, of the wake 

Syyfr, y, z) = 0 (2-3) 

Note that this surface oyy is considered twice, since a is a closed surface. In other 
words, the two sides of the wake are considered to be two independent surfaces 
having the same equation (but opposite outwardly-directed normals). 


*ln other words, the boundary conditions take info account the motion of the surface a, 
but are applied at the mean (steady) position. 
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Figure 2-1 . Aerodynamic Nomenclature for Potential Flow. 






As mentioned previously, the present formulation is based upon the Green’s 
function method. The infinite -space Green's function for the equation of potential is 
defined by 



= 6 (x-x*, y-y*, z-z^, t-tj 


(2-4) 


(where 6 is the Dirac delta function) with G = 0 at infinity. 

Multiplying the equation of the aerodynamic potential. Equation (1-2), by 
the Green's function G and subtracting Equation (2-4) multiplied by cp, yields 

(2-5) 
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dt 2 j 

CO 

\ 03 

/ 



Making use of the identities 


V ‘ (a?b) a Va * Vb + aV z b 


and 


dt 


db 

dt 


da 


db +g d 2 b 


dt dt 


dr 


Equation (2-5) reduces to 


V * (G?cp - cp^G) - — — 

2 dt 
*00 


f G & -cp — 
dt dt 


cp6 


(2-6) 


(2-7) 


(2-8) 


Multiplying Equation (2-8) by the domain function E, defined by Equation (2-1), 
integrating over the whole four-dimensional space-time and noting that, for any 
function g, (see Equation (2-1)) 


jjf £ 9 dV = J"J*j" g dV 


(2-9) 


-oo 


2 -4 



one obtains 


7 d * l 

v * (GV(p - <pVG) - — — (c ^ - tp — ) 
2 dt V dt dt/ 

-co V 

a 

L oo J 


dV 


Using the relations 


co 

= - 111! E » sdvdt 


( 2 - 10 ) 


-go 


oo 


J 


-oo 


Sg 

bt 


dt = g (co) - g (- 00 ) = 0 


(2-11) 


for g ( ± 00 ) = 0 and 


as 


JJJla dV = (JjJigS-Tdo =^g |VSy I” 1 do (2-12) 


V cr 


and similarly 

JjJ V - a dV =(jjj)a ' n da = (£j)° * v Sy|vSy ! ^ da (2-13) 

V a o 

where a is an arbitrary three-dimensional vector and n is the outwardly-directed* 
normal to the surface 


n = 


VSy 

!vs v ! 


(2-14) 


"Outwardly" is understood as "going from the body into the fluid", that is, from the 
region E = 0 into the region E = 1 . 
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Equation (2-10) reduces to 


j# 

-oo a 


VSy.(GVcp - cpVG) 


1 

_ L u — 

o oo _ 

2 dx 

a oo 


(g *£■ 

dG \ 

-cp 

1 *" 

dt /J 


oo 


- JjlJ* E cp 6 dV dt 

-oo 


(2-15) 


Finally, making use of the definition of the Dirac delta function yields 


oo 

E(p*)cp (p*/ t *)= J ^ 


-oo 


V Sy . (G^ cp - cpVG) 


± u (g 4 

2 dx l dt dt j 
a oo 4 


l-l 


|VSw| ‘ da dt 


(2-16) 


Equation (2-16) is the desired Green's Theorem for potential subsonic and supersonic 
small -perturbation aerodynamics. (Compare Equation (3-19) of Ref. 3.) 

2.2 SUBSONIC INTEGRAL EQUATION 

Equation (2-16) is valid for both subsonic and supersonic flow. In this sub- 
section the flow is assumed to be subsonic, i.e., M = U / a < 1 . In this case the 

' CO' CO 

Green's function is given by (Ref. 3) 

G = “ — — - — 6 (2-17) 

4 tt r . 6 
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where 


r 


P 




(y-y *) 2 


+ 



with 

p = >fl -M 2 

while 

6 0 =& (t-t* + 0) 

with 



fg + M (x-x*) 


(2-18) 


(2-19) 


(2-20) 


(2-21) 


2.2.1 SUBSONIC GREEN 'S THEOREM 


Combining Equations (2-16) and (2-17) one obtains 



4rr E(pJ cp(p*,t*)= - J ^IvSy.Vtp--^- ~~ ^ 

a L ^ SX * 


~6 |9S v r 1 d0dt 
9 


oo 


/# 

-oo a 


7S y .V 


<£, 3x ax\r p 


I’M 1 djdt 


OO 


-co a 


VSy . y6g 


U oo 3S V d6 0 


2 9x dt 

a oo 


<p — |vSy| ^dadt 


•P 


( 2 - 22 ) 
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Next note that 


36, 




at 


ve 




(2-23) 


Combining Equations (2-22) and (2-23) yields 


4tt E(P*) 9 (p*/ f *) = 


oo 


-j # 

-oo a 


■V’’’ ' 


IL, ss., j 

oo V dtp 

2 dx dt 

%o 


— 6 0 !vs y | 1 da dt 


oo 


J 0 

-oo a 


VSy-Vl 


U 2 3S., i 
oo V 3 

2 dx dx 

a oo 


96 0 IvSyl 1 da dt 


oo 


1 0 

-oo a 


3! 


V'-— j 1 + u ®f- 

a oo \ / 


36n , 

-2L — |VS.,| -1 dadt 

r at v 

p 

(2-24) 
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Performing the integration with respect to time, one obtains* 


4trE(p*) cp (p*, t*) = 



Equation (2- 25) is the desired subsonic Green's theorem. 

2.2.2 GENERALIZED PRANDTL-GLAUERT TRANSFORMATION 

In order to obtain a formally simpler expression for Equation (2-25), con- 
sider the generalized Prandtl-Glauert transformation and nondimensionalization 

$ = v Ay 

x = x/f3 i 
Y = y/JL 
Z = z/j l 

T = /i (2-26) 

where A is a reference length. 


In Equation (2-25), the symbol [ 


indicates evaluation at time = t* - 0. 


[ ] 9 h[ ] 
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Note that 


1 da 1 dy dz _ t dYdZ 


2 

jT dS 


r p IvSyl r p la^/Sxl R las^ysxl r |v o s v | 


(2-27) 


where v is the gradient operator in X,Y, Z variables, and E is the transformed a in the 
o 

space X, Y, Z, and R is given by 


where 


2 2 2,^/2 
r = [ (x-x*r + (v-Y*r + (z-z*r ] 


(2-28) 


Using Equations (2-26) and (2-27), Equation (2-25) reduces to 


4ttE(P*)§ (P*,T*) = 


-# CvvvTt 


dE 


R I V S 
I o V 


+ #Vv- 7 oWw 8 dz 


R 


I Vv I 


-# 


vs..* V 0 

o V o 




d T 


1 dE 


R V S..I 
o V 


(2-29) 


[ ]® h [ ] 


T =T* -0 


(2-30) 


with 


and 


© = U od Q/JL= [ M (X-XJ + R] M/p 


0 = [ M (X*-X) + R3 M/p 


( 2 - 01 ) 


(2-32) 
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Introducing the normal derivative 


— = N . V s — — V S>, • 7 
3N ° v s | o V o 

o V 1 


(2-33) 


Equation (2-29) may be rewritten simply as 



where 

Y = J!i- (2-35) 

3N 


indicates the component of the nondimensional velocity in the direction of the normal 
N to the surface E of the X, Y, Z space (not in the direction of the normal, n, to the 
surface a of the physical space) and is known from the boundary conditions (see 
Section 4). The relationship between Y and \|[ is given by 


Y = N .VS 


' ( 

BSy 

3$ 

+ — - 

3$ 

BSy 
+ I 

3$^ \ 

IVI' 

3X 

3X 

3Y 

BY 

3Z 

3Z/ 

-M 2 ) 


Bcp 


3cp 

+ 3S V 

3cp 

3x 

3x 

3y 

3y 

3z 

3z 


u 


CO 


L^i 5 <p.(n-M 2 n T) = 

IV v l * U 


co 


llfvi 

l 7 o S vl 


n 


3cp 

3x 


(2-36) 
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Equation (2-36) may be rewritten in a more interesting form by noting that 
i can be decomposed into two vectors; the first one, i = n^ n, normal to E and the 
second one, i^, tangent to E , in the (i, n) plane 


n n + f 
x t 


Note that 


where 


1 


n 


l-n 2 ) 



X 



0 

- n 

n 

■ = 

n n } 

1 

X 

y 


x y 1 

0 


n 


n n | 

1 


z 


X Z ) 


f=T t ./|T t | 


is a unit vector tangent to E and 


, t = l'tl 


n 2.2 2 2 2 2 
(l-n )+n n +n n 
' x ' x y x z 


1/2 


(2-37) 


(2-38) 


(2-39) 


Note also that 



= & 0 - mV ) 1 / 2 


(2-41) 
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Combining Equations (2-36) to (2-41) one obtains 


1 2 2 1-2 - 2 ^/ 2 - 
Y = — - (1 -M n ) Vcp'/n-Mn n n + (1 -n ) t 

■ t X I XX X 

u oo / 


1 2 2 2 2 - 2 2 - 

= (1-M n ) 7<p ' (1-M n )n - M n (1 — n ) t 

y X X XX 

1/9 /l- 2 \l/2 

■ — — — — (1-M 2 n 2 ) 7 t -( — ] M 2 n ^ (2-42) 

u co I '-mV/ ‘ 


where s t is the arc length along the direction of the unit tangent t . Note that if 
n = 1, then Y = B'I'/Uqo' whereas if 

Mn =0 (2-43) 

x 

then 

Y = -1— i]r (2-44) 

U co 

The approximation introduced in Equation (2-43) implies small-perturbation from free- 
stream velocity for M > 0 and is consistent with the use of the linearized differential 
equation for the perturbation velocity potential of Equation (1-2). 

It should be noted that, in Equation (2-34) the surface E is assumed to be 
fixed with respect to the frame of reference. However, the effect of the motion of 
the surface is retained in the boundary conditions, which are considered in Section 4.* 
Also note that (consistent with the hypothesis of small perturbation with respect to the 
steady-state configuration) the surface of the wake is assumed to be the one of the 
steady-state case. It may be noted that since E is fixed. Equation (2-34) represents a 
*This point is discussed further in Appendix C. 
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linear operator. Hence, flow can be decomposed into steady and unsteady contributions, 
which are decoupled and can be evaluated separately. As mentioned previously, the 
flow is assumed to be steady-state for T < 0. Therefore, the unsteady parts of the 
potential and the normal wash, $y and Yy, respectively, are identically equal to 
zero for T < 0: 

§y h 0, Yy h 0 (T < 0) (2-45) 

These conditions are introduced to permit the use of Laplace transform. This point is 
discussed further in Appendix B.4. 

2.2.3 CONTRIBUTION OF WAKE 

In order to understand the nature of the aerodynamic operator. Equation 
(2-34), it is convenient to isolate the contribution of the wake. 

Note that, as mentioned in Section 2.1, the surface E is composed of two 

branches. The first is the (closed) surface of the body. Eg. The second is the (open) 

surface of the wake, Tyj*:r\ote that this surface is considered twice since 2 is a closed 

surface. In other words upper and lower sides of the wake are considered to be two 

independent surfaces, having the same equation but opposite outwardly-directed 

normals, N and N ., respectively. This may be expressed as 
U it 


S = E B + E u + P- 46 ) 

N u = -Nj (2-47) 


*Note ; In inviscid flows the wake has zero thickness. However, the method may be 
used for finite -thickness wake with appropriate changes in the wake boundary condi- 
tions. For a fuselage with a truncated base, inappropriate modeling of the wake as a 
closed wake yields a nonunique solution to the integral equation. This phenomenon 
is analyzed in detail in Appendix E. 
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Combining Equations (2-34), (2-46) and (2-47) yields 




Note that, according to Equation (1-22), Ai|r =0. However, in general 
AY/ 0. For, using Equations (1-22) and (2-42), one obtains 

1/2 

M 2 n (2-49) 

x ds t 
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Note that- the small perturbation assumption (required for the linearization of 
Equation (1 -2))implies that n^ - 0. Therefore, AY - 0 and hence Equation (2-48) may 
be approximated with 


4nE(P*) § (P*,TJ = - (jjj) [ Y] 0 -L ds E 


0 


[$i 


© _b_ 

BN 


S5 

3T 


® 


J 3®_ 

R BN 


dZ„ 


♦ If 




j[A#] @ B 

( 1 ^ 


3 A? 

( 3N u' 

l R I 




0 


1 3 © 


R BN. 


( 2 - 52 ) 


where Eg is the closed surface of the body, while Eyy is the open (i.e., one side) sur- 
face of the wake, and A§ is the potential -discontinuity across the wake evaluated in 
the direction of the normal, N (i.e., M s } - § „). It should be noted that the 

value of A§ is notan additional unknown, since Equation (1-18) may be rewritten, in 
nondimensional form, as 


ti (P f T) = A$ (P JE , T -n) 


(2-53) 


where n is the nondimensional time necessary for a particle in the wake to travel 
(along a vortex line r within the steady flow) from the point Py E (origin of the vortex 
line at the trailing edge), to the point P. Small perturbation hypothesis is used for 
the steady-state flow, and hence, n is given by the nondimensional form of Equation 
0-19), i.e., 

n = P (X - X TE ) (2-54) 


Equation (2-52) is the basis for SOUSSA-P 1.1. 
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2.2.4 BASIC EQUATIONS USED FOR SOUSSA-P 1.1 

It may be noted that, as shown in Appendix A, Equation (2-52) is valid 
also on the surface £ if the function E, defined by Equation. (2-1), is generalized as 

E = 1 outside £g 

= 1/2 on £g (regular point) 

= 0 inside Eg (2-55) 

If the point P* h (X*, Y*, Z*) is on the surface £, Equation (2-52) is an integral 
equation (with differential -delay dependence upon time) relating the values of $ and 

— \ x * 

on the surface £ to the values of the normal derivative with Acp given by 

ST SN 

Equation (2-53). This is the equation used in SOUSSA-P 1.1. (On the other 

hand, if P* is not on £ , Equation (2-52) is not an integral equation but simply an 

integral representation of the potential $ at any point of the flow field in terms of the 

pi x 

values of $ and and the normal derivative on the surface £ .) 

ST SN 

2.3 SUPERSONIC INTEGRAL EQUATION 

Equation (2-16) was derived without assuming the flow to be subsonic; there- 
fore, Equation (2-16) is valid for supersonic flow as well. The supersonic Green's 
function is also defined by Equation (2-4) which may be solved to obtain (see Ref. 3). 

(2-56) 
(2-57) 


6— «- V + \ 

*"i 1 


where 


fp = < (x-x*) 2 " 0 ' 2 


(y-y*) 2 + (z-z*) 2 


1/2 
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with 


p ■ = y m 2 - 1 

while 


(2-58) 


with 


6 + = 8(t-t* -W 1 ) 
9" 


(2-59) 


0- = — L 


a oc P 


,2 


M (x* - x) ± r'p 


(2-60) 


and finally 


H (P) = 1 


x* - x > p 1 




= 0 


-X<P' 




■y*) 2 + (z- 


(2-61) 


Repeating the same procedure used to derive Equation (2-52) from Equation (2-16), 
using p' instead of p in Equation (2-26), i.e., 

X=x/p'j l Y = y/& Z =z/l T^u^t/je (2-62) 

one obtains 

4 tt E(P *)§ (P *,T*) = - P^[V]® + [Y-] 0 j "^-dE B 

+ Ul 0 + [§] @ 
e b\ 



(Equation continued on next page) 
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where E(P*) is given by Equation (2-55), 


JL =-N— + N— +N — 
9N C X 9X y SY Z dZ 


(2-64) 


is the conormal derivative. 


m _ 6$ 


y = 


dN 


is the conormal wash, whereas 


(2-65) 


► I — 

V — 


R' = I (X - X*) 2 - (Y - Y*) 2 - (Z - Z*) 2 


1/2 


( 2 - 66 ) 


and 


H = 1 for X* - X > (Y - Y*) + (Z - Z*) 


< 


2 2l 1 ^ 2 


= 0 for X, 


-x< [i:y-yj 2 + (z-z*) 2 V /2 


J 

] 


(2-67) 
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Furthermore 



indicates evaluation at time T = T* - 0 with 

0* = [M (X* -X) ± R'j M/p' 

and finally 

er = [M (X -X*) ± R'J M/p’ 


( 2 - 68 ) 


(2-69) 

(2-70) 


Equation (2-63) is the basis for the supersonic SOUSSA and is presented here for the 
sake of completeness and for future reference. The supersonic option is not available 
in SOUSSA-P 1.1. 
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SECTION 3 


NUMERICAL SOLUTION OF THE INTEGRAL EQUATION 

Equations (2-52) and (2-63) fully describe the problem of linearized unsteady 
subsonic and supersonic potential aerodynamics around complex configurations. In 
order to solve this problem, it is necessary, in general, to obtain a numerical approxi- 
mation for Equations (2-52) and (2-63). The numerical formulation used in SOUSSA-P 
is derived in this section. For the sake of completeness a general formulation for an 
arbitrary finite-element representation is presented first. Then the zeroth order 
formulation used in SOUSSA-P 1.1 is obtained. 

3.1 SUBSONIC FORMULATION 

Consider first the subsonic integral equation given by Equation (2-52). In 
order to discretize the space integral operator over the surface E it is convenient to 

D 

* 

use a finite-element representation for the normal wash Y and the potential $. 

3.1.1 SPACE DISCRETIZATION 

Using a general finite -element representation, it is possible to write 

(Ref. 8), 

H 

Y(P, T-0) = Y h (T-0 h )N h (P) 
h = 1 
H 

«(P, T-0) = £ $ h (T " 0 h )N h (F > (3-D 

h = 1 


Finite-element representation is meant here in a very broad sense: actually any inter- 
polation formula of the type given by Equations (3-1) and (3-2) is consistent with the 
formulation presented here. However, for arbitrarily complex configurations, only the 
finite -element interpolation (including splines over patches) is sufficiently general to 
be of interest here. 
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where Y^ (T -0^) and (T - 0^) are time dependent values of Y and $ at the point, 
P^, on Eg at the time T - 0^ (where is the disturbance-propagation time from 
P^ to P*); furthermore N^(P) are prescribed global shape functions, obtained by 
standard assembly of the element shape function (see for instance Ref. 16). The points 
P^ will be referred to as nodes*; H is the total number of nodes on the body. For 
simplicity the same shape functions are used for § and Y, although this is not essential 
to the method. 

Next consider the integration over the wake. In order to facilitate the use 
of Equation (2-53), it is convenient to divide the wake into strips defined by (steady- 
state) vortex lines emanating from the nodes on the trailing edge. The strips are then 
divided into elements with nodes along the vortex lines. The potential discontinuity 
can then be expressed as 

N 

A$(P,T-©)=£ M n (T~0 n )L n (P) (3-2) 

n = 1 

where N is the number of nodes on the wake, (T — © ) is the value of A$ at the 

n n 

nth node P^^ on the wake at time T - 0 (where 0 is the propagation time from 
n n n r r 

P^^ to P*), and (P) is the global shape function relative to the nth node of the 
wake. Note that according to Equation (2-53) 

- « < S) < T - n n> ^ 


*These nodes used to define the functions Y and § do not necessarily coincide with 
the points used to define the geometry of the surface Eg . 
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where m = m(n) identifies the trailing-edge node which is on the same vortex-line as the 

nth node P^\ Furthermore, n is the time necessary for the vortex -point to be con- 

vected from the trailing-edge node to the wake node . It may be worth 

° m(n) n ' 

noting that A§ = I , - $ . where h and h identify the trailing-edge nodes 

m n y u Jo 

(upper and lower sides, respectively) on the body corresponding to the mth node on the 
trailing edge. Therefore, it is possible to write 


«m(n) = E S nh*h 

h = 1 


(3-4) 


where S h = 1 (S nh = -1), if h identifies the upper-side (lower-side) node P h on the 
body corresponding to the nth node P on the wake (i.e.,Pg coincides with the node 

^m(n) on th ® trai, ' ng edge )r and S n h = 0 ot,ierw ' se * Thus ' 

^nh “ if ^ = ^m(n) ' s ° n u PP er s ' de 

= -1 if P^ = is on the lower side of Sg 


otherwise. 


(3-5) 


Combining Equations (2-52), (3—1) and (3-2) one obtains 
2E(P*)$ (P*,T) = Z B h Y h (T-® h ) 


+ ^ G h # h( T " e h5 + ^ D h i h (T “ 0 h ) 


+ S F n ^§ n (T-0 n )+ Z)G n Af n (T-0 n ) 


*For notational simplicity T* is replaced with T in this whole section. 
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where 



and according to Equations (3-4) and (3-5) 


(T-eJ = £s nh i h (T-0 n -n n ) 


n ' n 


n n' 


(3-8) 


Next consider, in particular, that P* coincides with the node { on Eg (P* = Pj). In this 
case E = 1/2 and, using Equation (3-8), Equation (3-6) reduces to 

5 i (T) “ B jh h ' ®jh' + ? C i h * h " V 

+ ^ D |h * h (T " V F in s nh * h (T ‘ 8 jn ' 11 n> 


Sj G. S,§,(T-0. - n ) 

n h |n nh h |n n 


(3-9) 


where 


B jh' C jh' D jh' F ]n' G jn' 0 jh “ [ B h' C h' D h' F n' G n' ®h 


P* = P. 

I 

(3-10) 
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3.1.2 


LAPLACE-DOMAIN ANALYSIS 


Equation (3-9) indicates the nature of the aerodynamic operator relating 
potential and normal-wash as obtained by using f ini te -element representation to dis- 
cretize the spatial problem. The operator is a linear differential -delay operator to 
which the methods of operational calculus can be applied. In this section $ and Y 
will indicate the Laplace transform of the unsteady part of the potential and normal 
wash, respectively, with initial conditions obtained from Equation (2-45).* However, 
before considering the Laplace transform of Equation (3-9), it is convenient to make 
some remarks about the contribution of the wake. It may be noted that, according to 
Equation (2-45), §y is identically equal to zero for T < 0; therefore, according to 
Equation (3-4), 


(«u> n = 0 a < n n ) 


(3-11) 


Hence, if the analysis is limited to T < T , only the contribution of 

elements with II < T need be considered, because II = T represents the 

n max n max r 

full extent of the physical wake. The elements with n > T would contribute to 
r 7 n max 

the transfer function and thus to the transform of $ but not to the final solution in the 

time domain for T < T . The advantage of deleting these latter elements is not 

max 

only that less computational time is used (since fewer elements are employed) but also 

-p n n 

that the issue of convergence connected with the infinite wake (factors e with 

Real (p) < 0 and n -»■ oo) need not be addressed.** 


*See Appendix B.4. 

**This problem is analyzed further in Appendix B.2. 
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Next, taking the Laplace transform of Equation (3-9) yields* 


L Y iK 


h - L- jhj "h 


(3-12) 


where § ^ and Y^ are the Laplace transforms of the unsteady parts of and Y^ whereas 

-p®. L 

Y ih = 6 |h " (C jh +P D |h ) 6 ' 


-p( 0 !n + n ) 


■E( F in + PV® " s 


nh 


(3-13) 


and 


Z .i — B .1 
|h |h 


-p® 


ih 


(3-14) 


In Equations (3-13) and (3-14), p is the complex reduced frequency (non-dimensional 
Laplace parameter), given by p = (y + I) k, where k = (Ujfc/U^ is the reduced fre- 
quency and y = Real (p)/ lmag(p). Note that for simple harmonic motion, p = ik 
(see also Appendix B). In particular, steady-state is obtained with p = 0. 

Equation (3-12) is the desired numerical approximation of the subsonic 
integral equation relating (in the Laplace domain) the transformed vector of the 
potential ^* e transformed vector of the normal wash |y^ j • It is the form 

that is actually solved in the SOUSSA-P 1.1 program. 

3.2 SUPERSONIC FORMULATION 

Using Equation (3-1) and following the same procedure used for the sub- 
sonic case, one obtains the supersonic time-domain aerodynamic operator 


*The same results are obtained for oscillatory flow (see Appendix B.3). For a dis- 
cussion on the initial conditions see Appendix B.4. 
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where 


2 E(P*) = (P*, T) = J2\ K (T " 0 h ) + V T “ 0 P I 
h 


+ E c h k(T-®h) + § h (T " 0 P 

h 


+ E[ D hM T -$ + Vh< T -^' 

h 

+ Vf' T&§ (T - e 4 ) + a$ (T - 0 ") 1 

Z-# n [_ n n n ' rr J 
n 

+ \ "* [g &f (T - 0 ) + G A§ (T - 0~)"j 
Z-/ [_ n n ' rr n n ' rr j 

n 



(3-15) 


(3-16) 
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with 


R' 


(X - X*) 2 - (Y - Y*) 2 - (Z - Z*) 2 


1/2 


(3-17) 


In particular, if P* coincides with the node j , using Equation (3-8), Equation 
(3-15) reduces to 

*i (T) =E B ih [v T - 8 iJ> + V T - 8 ji,>] 

h 

4 E c ih [ $ h (T " @ ih )+$ h (T "V] 

h 

+ E [ D {h 1 h (T ~ ®ih ) + D jh ? h (T ■ v] 

h 

+ EE F in S nh [ 5 h (T " 8 in - n n> + »h ( T - % ‘ «„>] 
h n 


where 


EE [°h S „h } h < T - e in - V + G i„ S nh ^ < T * % ‘ "„>] 

(3-18) 


h n 


' • i + 1 + + \ 

B.i , C.. , D. u , F. , G. _ © IT 1 
L |h |h' \n |n |n' jh/ 


B h' C h' D h' F n' G n' ®h 


P = P 
* * * • 


(3-19) 


Finally, taking the Laplace transform of Equation (3-18) one obtains* 


*The same results are obtained for oscillatory flows (see Appendix B.3). For a 
discussion on the initial conditions see Appendix B.4. 
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where 




Y ih = s jh - < c T h + p D f h > ^ pS ' h -< c ih +p V e 




- - -p 0 fh 


- V (F. + p Gt ) e 
“ P * I" 


-p(ot + n ) 

i n rv 




-p(e: + n ) 

— . * ' m n' 


-E(f! + p G: ) e * n 

“ I" l n 


nh 


nh 


~ i "P©;l 

Z.. = B. u (e l h 

|h |h 


+ e 


-PSjh) 


(3-21) 


Equation (3-20) is the desired numerical approximation of the supersonic integral 
equation relating (in the Laplace domain) the transformed vector of the potential 
to the transformed vector of the conormal wash { Y, ) . Note that, for 






supersonic trailing edges, F.^ s G. n = G.^ s 0, i.e., the contribution of the wake 
is identical!/ equal to zero. Note also that even if the trailing -edges are not completely 
supersonic (i.e., if a portion of the aircraft is contained in the Mach aftcone of any 
trailing-edge point), the wake can be truncated at finite distance from the aircraft, 
without any effect on the solution. 


3.3 


FORMULATION FOR SOUSSA-P 1.1 


This section deals with the subsonic formulation used in SOUSSA-P 1.1. 

This formulation is a particular case of the general formulation presented above, and 
is briefly illustrated here. Divide the surface of the aircraft Eg into small 
elements, E^ . Consider the shape -function, N^, equal to one inside E^ and equal to 
zero outside E^, i.e,. 


N h (P) = 1 


if 


p e Ei 


= o 


otherwise 


(3-22) 


3-9 



A point located on the element and identified as the center of the element will be 
designated as the point at which $ and Y are evaluated. Equation (3-1) may thus be 
interpreted as saying that, within the element , the normal wash and the potential 
are approximated with the values Y^ and at the center, , of the element, . 
(Note that the shape functions given by Equation (3-22) may be called zeroth-order 
shape functions. Therefore, the formulation presently used in SOUSSA-P 1.1 may be 
called zeroth-order finite-element formulation.) 

Next, note that using Equation (3-22), Equation (3-7) yields for instance 


B h 



C 


h 


1 

2 rr 



(3-23) 


(3-24) 


D h 


2 IT 



A 



(3-25) 


If is a quadrilateral element, then may be approximated with a 
hyperboloidal element (see Ref. 7) 

P = POO + ^10 ^ + PqI ^ + ^11 ^ 


1 1 

= E E F mn ^ 1 ' n H * W 1; -I S US 1) P’ 2 *) 

n = 0 m = 0 
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where E and T] are local coordinates with origin at the element center, P„ rt , and P 

' 00 mn 

are obtained in terms of the locations of the four corner points as 



(3-27) 


Note that quadrilateral hyperboloidal elements can be combined to yield a closed sur- 
face. Using these elements, the coefficients and can be evaluated analytically 
as (see Ref. 17) 


q, = I D 0,1) - l D (1,-l) - l D (-l,l) + l D (-l,-l> (3-28) 

Bfc = l s (l,l) - l s (l,-l) - l s (-l,l) + l s (-l,-l) (3-29) 

with [using -tt/2 < tan * ( )< rr/2 ] 

r 

I d (§,T1) = (1/2tt) tan p 1 (R x * Rx Aj/l R 1 R • Aj x A^) (3-30) 
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and 


lc(S,T|) = ~ — 

5 2 IT 


- R x A, ' N — - — sinh ^ 
1 c 


1 



where 


+ R x 


1 


sinh 


-1 




+ (R . N) tan 
' c p 


R x A-j • R x A ? 
| R| R • Aj xA 2 


R(SfT i) = p - p* = p 00 + P T0 § + p 01 h + P, } u\ - p* 

^(5,11) = 3R/3? = P 10 + T] P n 
A 2 (5,T1) = BR/BT1 = P Q1 + § P n 

N(5, Tl) = Aj X VIA, x \ I 


and the subscript c indicates evaluation at the element center (? = T| = 0). 


(3-31) 


(3-32) 

(3-33) 

(3-34) 

(3-35) 


If the values of ft * N at the four corner points of the element do not have 
the same sign, then the use of the principal value of tan ^ is incorrect (Ref. 18). 

In this case the tan ^ term Is evaluated by replacing the hyperboloida! element with 
two triangular elements: this is legitimate in view of the fact that the integrand in 
Equation (3-24) is the solid angle (see Ref. 3) which depends only upon the perimeter 
of the element, not the actual surface shape. A more complete analysis of the problem 
is given in Ref, 18. 
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Next consider the coefficient which is evaluated as 


1 

2 TT 


JJ i 


30 


dS, 


R 3N 


M 

2tt|3 


JM 


3R 


R \3N 


-MN X dz h 


M 

2np 


JH 

S h R ‘ 


dSi 


— dSi + M Ny II — Uil 

2 3N h 2tt[3 *h JJ R h 

S h 


'><!. JI » 


“ ‘ MN x h B h> M /P 


(3-36) 


where R^ and N^ are the values of R and = N • i at the center, P^, of the element 
^ h 

V 


Similarly the wake is divided into elements, E , and the shape -functions, 
L n (P), are defined as 


L (P) = 1 if P € s' 
n ' n 


= 0 otherwise 


(3-3 7 ) 


Hence the coefficients F and G are defined as 

n n 


F 

n 



Mds' 

SN U " 



A 



(3-38) 
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The definition of F n is identical with the one of (with E^ replaced by E 1 ). There- 
fore, the evaluation of F is identical to the one of C L . For G , assuminq that the 
wake has small inclination with respect to the X-axis (i.e., N^<< 1), the expression 
is 


G = R F M/p (3-39) 

n n n ' v ' 

where R is the value of R at the center of element E* . 
n n 

Finally 

®h = [ M( Y x *> + l f h - f *l ] m /p (3-<o) 

is the propagation time from the center Pj of the element E^ to the control point P* . 
Similarly 

®n = [ M < X n CW) - X .> + l F „ W " P.l] M /fi (3-41) 

and 

„n = , [xW - X h(n) ] (3-42) 

where h(n) indicates the center of the trailing -edge element from which the wake strip 

«~(W) v 

containing the point ' emanates. In the User's Manual (Reference 15), © n indicates 

e = e + n (3-43) 

n n n ' ' 
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SECTION 4 


NORMAL WASH 

In this section the relationship between the normal wash ¥ (or the conormal 
wash ¥') and the generalized coordinates, q^, describing the deformation of the air- 
craft, is presented. It should be noted that in this report the motion of the aircraft is 
assumed to consist of very small (infinitesimal) perturbations around a steady-state con- 
figuration which itself generates a small perturbation from free stream. In particular, 
in defining the integral equation, the surface. Eg, of the body was assumed to be time 
independent (see Equation (1-10) with e = 0). However, as mentioned in Subsection 
2,1, the unsteady contribution is retained in the boundary conditions which are con- 
sidered in this Section, This point is analyzed in detail in Appendix C where it is 
shown that in order to take into account the motion of the surface, ¥ in Equation (2-34) 
is given by Equation (C, 18), i.e.. 



where the superscript (t) indicates evaluation with time-dependent surface (e ^ 0 in 
Equation (1-10)). The vectors a & are evaluated from the unperturbed surface (e = 0 in 
Equation (1-10)). The subscript M on ^ indicates that the modified expression for t 
(including the motion of the surface) is used here (see Appendix C). 

4.1 PHYSICAL BOUNDARY CONDITION 

The boundary condition on the body is obtained by imposing the condition 
that a fluid particle which is on the point p of the surface of the body at time t will 
remain on the surface, i.e., will be on the point p + Ap (also on the surface of the 
body) at time t + A t . If the surface of the body is described by 

S (p, t) = 0 (4-2) 


4 - 1 



then it is also true that 


S (p + Ap, t + At) = 0 


(4-3) 


with Ap - V pAt, where Vp is the velocity of the fluid, Vp = V 0 . Taking the Taylor 
series of Equation (4-3), and using Equation (4-2), one obtains 


— +VS • V c I At + O (At 2 ) = 0 
dt *' 1 


(4-4) 


or, taking the limit as At goes to zero. 


DS 

Dt 


= 0 


(4-5) 


where, using Equations (1-1) and (1-3) 


— = — + Vp • V = — + y 0 . v 

Dt at r at 


a A a A d , 

— + u — +V<p. 7 = — + V 9 ' v 
at 00 ax dt 


(4-6) 


Note that 7S can be assumed to be directed as the outwardly directed nor- 
mal n (eventually Equation (4-2) must be multiplied by -1) so that 


n’M = «- 

|VS| 


(4-7) 


The superscript (t) emphasizes that the normal is evaluated at t since n indicates the 
normal in the steady-state configuration. 

Thus, the boundary condition. Equation (4-5), may be written as 

<!£ + ■ 1L = o 


|vs| 


|vs| dt 


(4-8) 



or 


^M 



1 

|7S| 




oo 



(4-9) 


4.2 RELATIONSHIP BETWEEN NORMAL WASH AND DISPLACEMENT 


Next note that while Equation (4-9) is the classical form for the boundary 
conditions, here it is convenient to express the boundary condition in terms of the dis- 
placement, u. Therefore, the surface of the body will be assumed to be given not by 
Equation (4-2) but rather as 


X - x(t) - p if) + t) (4-10) 

where p gives the reference configuration rigidly connected to the i, j, k-space while 
u is a small steady or unsteady displacement relative to this configuration. In Equation 
(4-10), E? y (a = 1,2) are convected curvilinear coordinates, moving with the surface 
of the body. The normal n^ is then given by 


n-<‘> = 


=-<’> X of 

if x of I 


(4-11) 


with 


a 


9x _ 5p 


+ 




+ 


u 

a 


(4-12) 


where u i^du/c^. 
a 

In order to rewrite Equation (4-9) in a form compatible with Equation (4-10), 
Introduce the velocity of the surface of the body 
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(4-13) 


V B = 


Sx _ du 


at at 

Then the component of v R along the normal n^ is given by* 


; B ' n ' (,) = 


Note also that the x-component of the normal n^ is given by 

n (t) = as L s | 

x 3x / 

Therefore Equation (4-9) may be rewritten as 

J + v P 


M 


oo 


v B ]*nW 


(4-14) 


(4-15) 


(4-16) 


Finally combining Equation (4-1), (4-11) and (4-16) one obtains 


Y = 


U 


U oo 5 + V B 


oo 


a-f) x S<*> 


l a i x a 2 | 


(4-17) 


where 


Note that using Equation (4-12) one obtains 


a-('> x of 


I°1 x a 2l 


n + An 


(4-18) 


An - (u.j x a 2 + a j x u 2 + Uj x u 2 j/ |a^ x a 2 | 


(4-19) 


Equation (4-14) may be obtained from Equation (4'r5) written for a_point moving with 
the surface, i.e., having velocity v D . This yields /\V S B | + v,, • n = 0 in agree- 

b ~5T/ b b 

ment with Equation (4-14). Note that Equation (4-13) is valid only for a wind-axis 
formulation. The extension to bodyAaxis formulation is given in Appendix D. 
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This yields 


Y = [ -7 + — L u J • (n + An) = Y $ + Y y 

^oo 


(4-20) 


where the steady part of the normal wash is given by 


Y s = -i . n 


(4-21) 


and the unsteady part is given by 


U 


L U ao 


u • (n + An) - i • An 


(4-22) 


4.3 LINEAR RELATIONSHIP BETWEEN NORMAL WASH 
AND GENERALIZED COORDINATES 


Let the displacement u be expressed as 
M 

5K".t) = %,(*> “■»<*"> < 4 - 23 > 

m = l 


where q m are the nondimensional generalized (Lagrangian) coordinates and M m (5°) 
are prescribed nondimensional modes (ordinarily the natural undamped modes of vibra- 
tion of the aircraft). In SOUSSA-P 1 .1 the motion is assumed to consist of small 
perturbations around a steady-state motion. Therefore q m are assumed to be small. 
Then, to first order one obtains 


A n 



•m 



(4-24) 


I _ _ BM \ 

An m = -M x a 2 + a l x /l°l x a 2l ( 4 " 25 ) 

\ SI 1 ss 2 / 
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and, again neglecting higher order terms. 


M 


M 


U 


U oo 01 = 1 


2 M m *nq m - E i'An 


m = ] 


q 

m m 


Next, taking the Laplace transform of Equation (4-26), one obtains 

M 


’’h = £ 

m = 1 


sM • n - i • An 
U m m 

co 


■m 


P=P h 


(where s is the Laplace parameter) or 


:0) 

"hm 


'm 


with 


c0) 


nm 



pM * n - i • An 
r m m 


J P = Ph 


where 


P = sA / U co 

is the nondimensional Laplace parameter (complex reduced frequency). 


wash matr 


Equation (4-29) is used in SOUSSA-P 1.1 for the evaluation of the 
"< [" E hm] • 


‘It is assumed q n (0) = 0 (see Appendix B.4). 


(4-26) 


(4-27) 


(4-28) 


(4-29) 


(4-30) 


normal - 
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SECTION 5 


PRESSURE COEFFICIENTS AND GENERALIZED 
AERODYNAMIC FORCES 

In this section the relationship between the pressure coefficient and the 
potential is considered first. For the sake of completeness the formulation for arbitrary 
finite -element representation is included here. The modifications necessary for the 
zeroth order formulation used in SOUSSA-P 1.1 are then outlined. These include a 
special definition of the shape functions (making use of the averaging scheme first 
introduced in Ref. 9) and the special evaluation of the pressure coefficient for trailing- 
edge elements (for both subsonic and supersonic flow). Finally the relationship between 
generalized aerodynamic forces and pressure coefficient is presented for an arbitrary 
as well as zeroth-order finite-element representation. 


5.1 PRESSURE COEFFICIENT 

The pressure coefficient is evaluated from the linearized Bernoulli theorem as 


- ^ dcp _ _ 2 l 


U 


2 dt 


U 


co 


CO 


Sf y 5$ 
^ °° ^ 


(5-1) 


In order to evaluate — , it is convenient to use the following procedure. Note that 
dx 


3$ 

— 3 V? 
dx 


(5-2) 


and that y $ may be evaluated as 


y $ = i±- a 1 + iL a 2 + n 
d ? 1 d ? 2 Sn 


(5-3) 


5 - 1 



where 


with 



The covariant components 


of the metric tensor are given by 


°12 °2l) 


(5-4) 


(5-5) 


«0 


a ♦ 
a 


a 


0 


where the base vectors a^ are given by 


(5-6) 


a 


dp 

df 


(5-7) 


Therefore, 


= -21 



3$ 

at 


1 

d5 


-2 


a 


nii\ 
h 2 9n / 


(5-8) 


Next note that, as with Equation (3-1) 

H 

$ (f , T) = E «u(T) N. (f) 

h = l h h 


(5-9) 


(where (§ a ) does not necessarily coincide with N^). Hence, the tangential deriva- 
tives are evaluated as 



H 


h = 1 


3N. 


a? 


. 0 / 


(5-10) 
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Thus neglecting the term n — and combining Equations (5-8) and (5-10) one obtains 

x bn 


C <p k , t) = -2 -5- — -2 E T • 

P 4x, h 


-1 aN h , -2 SN h 
a + a 




l h (5-11) 


55 h~n 

P-P k 


or in the frequency -domain 

C S.) =-2pi k -2 Er- 


5Ni bN, 

-1 . -2 

a , + a ~ 




1 


sr 


(5-12) 


P=P k 







c , 

p,k 


1 

-C 

LL1 

l 


? h 


with 


(5-13) 
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a’ 1 


bN 


^ + a 2 





P=p" k 


FORMULATION FOR SOUSSA-P 1.1 


(5-14) 


Next it should be noted that the formulation for the evaluation of the pres- 
sure as described in Section 5.1 needs modification in the case of zeroth-order form- 
ulation, for Equation (5-10) cannot be used in connection with Equation (3-22) unless 
it is interpreted in terms of the theory of distribution (or generalized functions). In 
order to avoid this problem, the following procedure is used in SOUSSA-P 1.1; from 
the values, of the potential at the centers P^ of the elements, E^, the values, 

I _ 

$ 1 ^, of the potential at the corners, P£, of the elements, are evaluated by averaging. 
In other words 
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(5-15) 


( ’k 5 = IA kh> [ V 


where [ A^J is a weighted averaging matrix defined as 



if 


p k 6 Hi 


(5-16) 


(i.e. , if the P£ is one of the corner points of the element £^) and 

Akh = 0 otherwise. (5-17) 

In Equation (5-16), the weights are proportional to a typical length of 
the element (which is assumed to be equal to the square root of the area of the element 
E^). Having evaluated the values of $ at the corner points, the potential is expressed 
as 


* = E ^ (5 1 ,? 2 ) (5-18) 

k 

* 

where are first-order global shape functions obtained by assembling local shape 
functions of the type 


N 


(E) 

k 


1 

4 *k\ 


(5 + S k ) 01 + \) 


(5-19) 


where §k = ±1 and = ±1 are the locations of the corners of the element £^ 
(5 and T| are the coordinates over the element; see Equation (3-26)). 


Thus, Equation (5-10) is replaced by 



an' 


■ EE 


3N, 




a? 


k if 




(5-20) 
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Equation (5-20) is formally equal to Equation (5-10) with 


3N h 



E 





(5-21) 


However, note that the averaging scheme cannot be used to evaluate the 
value of the potential difference between upper and lower sides at the trailing-edge 
nodes. Therefore, a different method has to be used to evaluate the pressure at the 
centers of the elements adjacent to the trailing-edge. 

Consider the subsonic case first. In this case the trailing edge evaluation 
is based upon the Kutta condition, i.e., that the pressure discontinuity goes to zero 
like the square root of the distance from the trailing edge. Therefore, near the trailing 
edge. 


- k, 


+ k °/ 


Xy£ “ X 


C p x " k A " k D y x TE " x 


(5-22) 


Noe that the average value of the potential (cp^ + cp^)/2, at the trailing 
edge can be evaluated correctly using the averaging scheme. Hence, indicating with 
Cp the values obtained with the averaging scheme, k^may be obtained from the values 
of C* at p. and p. (assuming that x. - x. , see Figure 5-1), as 

p ij i 2 i 1 i 2 


C + C = constant = 2k. - C' + C' 

P a A P; P; 


(5-23) 


1 


Note that, according to Equation (5-22), k^ may be evaluated from the 

values of C' at p. and p. (assuming that x. - x. , see Figure 5-1), as 
P '3 '4 '3 ‘4 
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ELEMENT BOUNDARY (TYPICAL) 


Oi 

i 

O' 



Figure 5-1 . Geometry for Trailing Edge Evaluation 


WAKE 


“ Cp ) ( X TE “ x ) ^ = constant -2k 

r u r £ 


D 


«(C* - C ) 
P: P; ' 

'3 '4 


x. + x. 
'3 '4 


,-1/2 


TE 


(5-24) 


Equations (5-23) and (5-24) yield 


k n = I C' - C' 

D o 1 P; P; 

'3 '4/ 


X TE " 


x. + x. 
'3 '4 


1/2 


= - O + C' ^ 
2 P: P; 

1 \ 'l l 2l 


(5-25) 


Equations (5-24) and (5-25) may be used to evaluate C and C (assum- 

P ii P io 


ing that x. - x. ) as 
'1 '2 


C = 1 (c: +■ cl + pc' - PC' ) 

Pi, 2 Pi, Pi 2 Pi 3 Pi 4 


(5-26) 


c =!(£' + £■- *p£' + }c ) 
Pi 2 2 Pi, Pi 2 Pi 3 \ 


(5-27) 


where 


2x 


TE 


P = 


x. ” x. 
'l '2 


2x 


TE 


x. “ x. 

'3 '4 J 


1/2 


(5-28) 
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Note that Equation (5-26) together with 


C = O (E. not on trailing edge) 

Pj Pj 1 


(5-29) 


may be rewritten as 


hi -Klls, 


(5-30) 


(TE) 

where the matrix E:j ' is obtained from the unit matrix (see Equation 5-29) by replacing 

the row corresponding to the upper-side elements E. which are in contact with the 

'l 

trailing edge (see Equation (5-26)), i.e. f for i = ij, with 


E<T E) = 1/2 for j = i, 

= 1/2 for j = i 2 

= p /2 for i = '3 

= -p /2 for { = i 4 

= 0 otherwise 


Similarly for the lower-side elements E. , i.e., for i = L , 

'2 Z 

E<| E > = 1/2 for i = 1, 

= 1/2 for j = i 2 

= -p /2 for i = '3 

= p /2 for j = i 4 

= 0 otherwise 


( 5 - 31 ) 


(5-32) 


5-8 



This is easily accomplished in the program by introducing an index function 


'h “ 'Kutta < h ' m) 


(5-33) 


(where m ranges over all the trailing edge segments and h = 1, 2, 3, 4) which gives 
the element numbers ij, ^ ig and i^ relative to the m^ 1 trailing edge segment. 

The above results may be combined by noting that 


hf ] - k>] k>] 


(5-34) 


gives the proper values for C at all the element centers as 

P 


•4 ■[ 


; (3) 

fc kh 


(5-35) 


5.3 


GENERALIZED AERODYNAMIC FORCES 


Next consider the generalized aerodynamic forces given by Equation 
(1-13). In order to obtain the relationship between e^ and C^ ^ set, in analogy to 
Equation (3-1) 


H 


c p = E c P .h N iX> 

h =0 

Combining Equations (1-13) and (5-36) yields 


(5-36) 


e = 
n 


-1 


S? a 




C i N, n ’ M da 
p,h h | n 


(5-37) 
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Equal-ion (5-3 7) may be rewritten as 


<•„! ' [E S' lc P ,h ! 


where 


= — (ff)Ni n • M da 
nh o n n 


je 2 a 


(5-38) 


(5-39) 


In SOUSSA-P 1.1, in approximating Equation (1-13) and hence Equation 
(5-39), the pressure is assumed to be constant within each element, i.e., the pressure 
mode shapes are given in Equation (3-22). (This is consistent with the numerical approxi- 
mation of the integral equation.) Combining Equation (3-22) and (5-37) one obtains 



(5-40) 


or, approximating the integrand with its value at the center. 



( 5 ! x a 2 



(5-41) 


Note that a are the base vectors in the physical space, not in the Prandtl-Glauert 
space. It should also be noted that the formulation presented in this subsection is 
valid for both subsonic and supersonic flows. 
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SECTION 6 


APPLICATIONS 


In this Section the formulation developed in the preceding Sections is 
applied to two specific problems: aerodynamic analysis (with emphasis on flutter 
analysis) and aerodynamic design. 


AERODYNAMIC ANALYSIS 


For the sake of clarity matrix notations will be used throughout this 
Section. For instance, Equation (4-28) may be written as 


l = h q 


i=R} < 6 - 2 

9 = (U ( 4-3 

il = [ i hm ] (4- 4 

Similarly Equations (3-12), (3-20), (5-35) and (5-38) may be written as 
Y § = Z Y (6-5 


Sp = -3 »- 


e = E . C 
-4 -p 


6 - 1 



where r is given by Equation (6-2), $ - (® n 3 / C - [C^ and e (e^}, whereas 




(6-8) 

(6-9) 

(6-10) 


Combining Equations (6-1), (6-5), (6-6) and (6-7) one obtains 


i = ! 2 i 2 

£p = -321 2 

e = E q 


-21 = -2 -1 
-321 = ^3 -2 -1 
I = -4 -3 -2 ^1 


( 6 - 11 ) 

( 6 - 12 ) 

(6-13) 

(6-14) 

(6-15) 

(6-16) 


where 


E 2 = y" 1 z 


(6-17) 


6.1.1. PHYSICAL MEANING OF E 2] , E^ AND E 

Note that the nth columns of the matrices E^, £321 an< ^ ^ 9' ve ^* e vec t° rs 
of the values of the potential, pressure coefficients, and generalized forces, respectively, 
due to the nth boundary-condition mode (nth column of Ej). This is apparent from 
Equations (6-11), (6-12) and (6-13) by setting 
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q =1 m = n 

= 0 m / n (6-18) 

Finally, note that the matrix E given by Equation (6-12) is the matrix of the 
generalized aerodynamic forces used, for instance, in flutter analysis. (This is apparent 
from Equation (6-13).) Other applications for the matrix | include evaluation of 
aerodynamic coefficients, stability derivatives and forces due to turbulence. 


6.2 AERODYNAMIC DESIGN* 

Consider the aerodynamic-design problem (also known as the inverse problem) 
of determining the aircraft surface shape which will generate a desired (prescribed) 
steady-state (p = 0) pressure distribution. 

The problem may be solved using the procedure presented here. Let p^ 
and p^ indicate an initial-estimate shape and the revised shape, respectively. Let 


the displacement u = p^ - p^ be expressed 


as 


N 




1 


(6-19) 


The difference in pressure between initial and new shape is given by 



( 6 - 20 ) 


where (p = 0, for steady state) 


E 




( 6 - 21 ) 


If the pressure is prescribed at N points, solving Equation (6-20) one obtains the 
values q n and hence the displacement u and from this the new shape. If necessary 
the procedure might be repeated until the iteration converges. (At any iteration a 


*This option is not available in the version 1.1 of the program SOUSSA-P. 
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new geometry is considered; this implies re-evaluating the matrices Ej, E 2 and Eg, 
which depend upon the geometry.) 

Note, however, that the above procedure is likely to yield wiggly shapes. 
Hence it is convenient to use few mode shapes and satisfy Equation (6-20) in a least- 
square sense, i.e. * 


.T _T cT I / c I _ minimum 


6' -q’ E' I E q - 6 


( 6 - 22 ) 


i.e. 


A <pA A w 

e'e q = E 6 


(6-23) 


which yields 


-(Vi)-' V* 


(6-24) 


Equation (6-24) may be used to obtain the values of q from the desired 
(prescribed) change 6 , in pressure distribution. 


*Smoothing between iterations is an alternative procedure. 
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SECTION 7 


CONCLUDING REMARKS 


A general formulation for subsonic and supersonic, steady, oscillatory and 
unsteady aerodynamics for complex aircraft configuration has been presented. The 
formulation is used in the computer program SOUSSA-P 1.1. A few comments about 
the present formulation are presented in this Section. 


7.1 ADVANTAGES OF FORMULATION 


The formulation presented here offers a unified approach to solve steady, 
oscillatory and unsteady problems for both subsonic and supersonic flow. The formula- 
tions for subsonic and supersonic are very similar and highly "compatible" from a 
programming point of view. The formulation is such that very little human interven- 
tion is necessary; each surface element is treated the same way. 


The formulation is very modular. The matrices E^, Eg, Eg and E^ are ob- 
tained in a totally independent way. Note that the formulation is particularly advanta- 
geous for multiple frequency evaluation. For, note that, comparing Equations (4-29) 
and (4-30), it is possible to write 


E - E<°> + o I* 1 * 

ti " I] + P -1 


where E^ and E^ are independent of p. 


E = + D gO) 

b .3 £3 + P £3 


(7-1) 

Similarly (see Equations (5-14) and (5-34)) 

(7-2) 


where E^ and Eg^ are independent of p. Also (see Equation (5-39)), E^ is 
independent of p. Finally, 
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(7-3) 


~ ~ ~ 

§2 = Y Z 


where for subsonic flow Yj^ and Z.^ are given by Equation (3-13) where S.^, C.^, 

* 

Djh, e. h , Fj n , G.^, 0. n , n n and S^ are independent of p. Therefore, once the 
above p-independent arrays are evaluated, the evaluation for multiple frequency 
requires only the assembly of the matrices Ej, Y, Z and Eg and the evaluation of 


E - = E -4 -3 -21 


(7-4) 


where E 2 ^ is obtained by solving the system of algebraic equations 


Y E -2, ■ SI, 


(7-5) 


In addition, calculations for updated sets of modes (as necessary in optimal 
design where a different set of modes is used at each iteration, but geometry and Mach 
number are not changed) are greatly facilitated by the modularity of the formulation. 
For, only the matrices E ^ and E^ are mode -dependent, therefore, the matrices Y, 

Z and Eg need not be reevaluated at each iteration. 

The above features make the program SOUSSA-P 1.1 not only general and 
flexible, but also modular, simple to use and efficient, especially for calculation in- 
volving multiple frequencies or mode changes. 


7.2 WORK IN PROGRESS 

While the program SOUSSA-P is the most advanced program for unsteady 

aerodynamics for complex configurations, the Green's function method covers cases 

which are not included in the present version 1.1 of the program. In particular, 

the method can be used to solve nonlinear subsonic unsteady flow in the time domain. 
— ^ 

*Similar!y, in supersonic flow, Yj^ and Zj^ are given by Equation (3-21) where all 
the coefficients are frequency independent. 
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This work is now being completed. Preliminary results are presented in Reference 22. 
The formulation is now being extended to unsteady transonic flow (Reference 23). 

(The Green's function method is applicable to fully nonlinear unsteady transonic flow 
with moving shock waves (see Reference 3).) 

Additional work now underway includes a higher -order finite-element 
formulation, special-purpose elements (such as hinge-line elements), wake roll-up 
and jet engine flow modeling. Finally, preliminary work for the inclusion of the 
rotational part of the velocity flow field (for use in boundary-layer analysis) is also 
underway. 
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APPENDIX A 


THE VALUE OF THE FUNCTION E ON THE SURFACE 


Consider Equations (2-52) and (2-63) which give (for subsonic and supersonic 
flow, respectively) a representation of the potential §, anywhere in the fluid volume, 
in terms of the values of § and -J|q- on the surface of the body and 6? on the surface 
of the wake. The values of 3§/9N on the surface of the body are given by the boundary 
conditions but the values of § are not known. In order to solve the problem, it is thus 
necessary to obtain first the values of § on the surface. This can be accomplished by 
letting the point P* of the fluid volume V approach a point P q of the surface. Then 
Equation (2-52) yields 


where 



(A-l) 


(A -2) 


A - 1 



is the contribution of the wake and 


E* = 1 P* outside Eg 

E*= 0 P* inside Eg (A -3) 

A similar expression is obtained from Equation (2-63) for supersonic flow. 

In this Appendix, It is shown that, in the limit. Equation (2-52) is still 
valid if the definition of the function E is generalized as follows: 

E # = 1 P* outside Eg 

E* = 1/2 P* on Eg (regular point) 

E* = 0 P* inside Eg (A-4) 

In order to simplify the discussion of the limit as P* -» P , the steady subsonic 
case is considered first. The results are then extended to the unsteady subsonic case. 
The supersonic case is considered last. 


A.l STEADY SUBSONIC FLOW 


Consider the steady incompressible flow first. By letting P*-» P , the 
integrands become singular in the neighborhood of P Q . Thus, it is convenient to 
separate the contribution of a small neighborhood of P , which will be indicated as 
E g . The neighborhood E g is a small circular element of the surface Eg with center 
P q and radius e . Thus, Equation (2-52) may be rewritten as: 


4ttE*$* = 



(A -5) 


A -2 



where I ^ is the contribution of the neighborhood of P , given by 


I = - 
e 



(A-6) 


A. 1.1 REGULAR POINT 

Consider first the case in which P q is a regular point of the surface Eg, 
i.e., a point in which there exists a unique tangent plane to the surface Eg . For 
notational simplicity, the analysis is performed by assuming that the origin is located 
at P q and that the Z-axis is normal to the tangent plane and directed from E = 0 to 
E = 1, i.e., parallel to the normal N. The point P* is on the Z-axis and therefore 
P* s (0, 0, Z A ). (The extension to an arbitrary element is obtained by replacing the 

A A A _ 

coordinates X, Y, and Z with a local coordinate system X, Y and Z with P q located 

A A A A 

atX=Y = Z = 0 and with Z normal to the tangent plane and directed from E* = 0 
to E*= 1.) Separating terms of order e and noting that X* = Y* = 0, Z = 0 and 

, Equation (A-6) reduces to 

Z =0 


SR 

bZ 


Z=0 


SR 

SZ* 



- Y 


II 


X 2 +Y 2 < e 2 


J X 2 + Y 2 + Z 2 


dXd Y 


- K II — - dXdY + O(e) (A-7) 

2 2 2 aZ * / y 2 + v 2 ":72 

X+Y < & * VX +Y +z* 

where the subscript o indicates evaluation at P Q . By using polar coordinates 

R = /x 2 + Y 2 

0 = tan ^ ^ (A-8) 
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one obtains 


I = - 2rr Y 
e 




/r 2 + z 2 


RdR 


- 2rr $ 


/ 




7r 2 + z 2 


RdR + O(e) 


Noting that 


1 3 


Z* 9^* 


R 2 + Z 2 2 


/ 


R 2 + Z 2 


I JL 

A A 

R 3R 


Equation (A-9) becomes 


| = -2uT 

e o 




n e 


R 2 + zl 


- 2 TT Z* 


J 1 


2 2 
IT + ZZ 


+ O(e) 


Finally, by letting P* go to P Q ,(that is, Z -* 0), one obtains 


lim I = lim 

p*-pl z*-o 


" 2tt¥ o / e 2 + Z 2 - jZ*| 


- 2 TT $ Zj 




2 , y2 

e + Z* 


|z.| 


+ O(e) 


- 2 tt Y e + 2 tt § sgn (Z*) 
o o 


+ O(e) 


± 2 it $o + O(e) 


(A-9) 



(A-10) 


(A-ll) 


(A-12) 
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where the upper (lower) sign holds for Z* > 0 (Z* < 0), that is, when P* is located 
outside (inside) the surface Eg ; correspondingly, the function E assumes the values 

1 (E*= 0). 


Finally, using this result in Equation (A— 1 ) one obtains, for P* on Eg 


(i .e. , 5^ $*)/ 


4tt E * — $ 


- S 




Y - dE 
R 


♦ // 


JJ 9 dE+l w + O(e) 
Eg-E e 3N\R/ 


(A— 1 3) 


Note that, in both cases (P* inside or outside Eg), 


- 


P* outside Eg 


0+1 

2 


= — P* inside E r 


(A-14) 


Furthermore, R is the distance between the dummy point, P, and the control 
point (on the surface Eg), P Q . Flence, by letting e go to zero. Equation (A-13) yields 


(fK +| w 

Z B E B 


(A -15) 


Note that for steady subsonic flow. Equation (2-25) yields 


4ttE*$ 


*=-# 


¥ i dEB+ # § ™ u dE B + ' w 


(A-l 6) 
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Note that Equation (A-16) must be used if P* is outside or inside the 
surface, v/hereas Equation (A-15) must be used if P q is on the surface. However, 
by comparing Equations (A-15) and (A-16), it is easily seen that Equation (A-16) is 
valid everywhere (outside, inside and on the surface Eg), if the convention is made 
that E is given by Equation (A-4). 


A. 1.2 COMMENTS 

It should be emphasized that in obtaining Equation (A-15) from Equation (A-13) 

the limit e-» 0 is performed with P* on the surface Eg. This implies that the contribution 

of E is now of order e . In order to clarify this point, consider the quantity 
€ 



and note that 


lim 

z*->o 





= -sgn(Z*) 


(A-1 8) 


whereas 


ii 


im 


e — >0 
o 


( lim | e ) 

Z.-0 e °! 


= li 


im 

e ->0 
o 


0=0 


(A-1 9) 
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The difference between these two limits is due to the fact that, in the limit (as Z*~»0), 

the integrand of l e behaves like a Dirac delta function and hence, its contribution 
e o 

for a domain which excludes the singular point is zero. 

It is apparent that the sequence of limits as indicated in Equation (A-18) is 
performed for Equation (A-12), whereas the one indicated in Equation (A— 1 9) is used for 
Equation (A-15). The above results may be restated by saying that the limit as P*~* P q of a 
doublet is a generalized function, i.e. 


lim 


-1 


-1 


p *-* p 0 aN \ 4ttR / 


BN 4ttR 


P*= P_ 


+ 1 8(P -P ) (A-20) 

2 ° 


where the upper (lower) sign holds for P q outside (inside) Eg. 


A.l .3 NON -REGULAR POINT 

The above results may be generalized to non-regular points (i.e., points of 
E S j in which there is not a unique tangent plane to Eg) if they are expressed in terms 
of the classical concept of solid angle. Note that by definition, the solid angle, dQ, is 


dQ 



dEg 

- cos a RN 

R^ 


R • N 



a 

BN 



(A— 21) 


where dEp is the projection of dEg into a plane normal to R (see Figure A-l). 
Using Equation (A— 21 ) one obtains 



0dQ = Q*/4 tt 
e b 


(A-22) 


A - 7 




where 


fi* = 0 P* inside Eg 

= 2tt P* on Eg (regular point) 

= 4rr P* outside Eg (A-23) 

Next, note that using Equation (A-22), Equation (A-16) may be written as 

4tt§* = - (ffi Y — dEg 

4 ? R 

S B 

(*-«*) — 

SN 

since, according to Equations (A-22) and (A-23) 

4ttE* " (“) dX B = 4itE * + Q * = 4 " (A-25) 

e b 

for P* outside and inside Eg. Note that the right-hand side of Equation (A-24) is 
continuous because the Dirac delta function in Equation (A-20) is now multiplied by 



ds B + *W 


(A-24) 


(*-**) Ip = p* = 0 (A-26) 

Therefore, Equation (A-24) is valid for all cases, i.e., for P* outside Eg, inside Eg 
as well as on Eg (for both regular and non-regular points). Finally, Equation (A-24) 
may be rewritten as Equation (A-16), if E* is defined as 



(A-27) 


A -9 



Note that the above Equation is a generalization of Equation (A-4) in that it is valid 
for non -regular points as well. 

In SOUSSA-P the control point is always located. at a regular point (center 
of the hyperboloidal element described in Subsection 3.3). Therefore, in the rest of 
this report it is always assumed E = 1/2 on Eg. 

A. 2 UNSTEADY SUBSONIC FLOW 

The results of the Subsection A. 1 .3 are immediately extended to unsteady 
subsonic case by noting that, using Equation (A-25), Equation (2-52) may be rewritten 
as 



(A -28) 


with lyy given by Equation (A -2). Equation (A -28) is valid for P* outside and inside 
Eg. In addition, there are no discontinuities as P* approaches Eg since the Dirac delta 
function due to the doublet distribution is multipled by zero. Therefore, the equation 
is still valid in the limit as P* approaches the surface Eg. This implies that Equation 
(2-52) is valid on the surface Eg as well if E* is defined by Equation (A -27), or in 
particular, by Equation (A -3). 
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SUPERSONIC FLOW 


A. 3 

Consider the steady supersonic case. For steady state Equation (2-63) re- 


duces to * 


4ttE*$* 



ds B 



where 1^ is given by 




(A -29) 


(A -30) 


A. 3.1 SUPERSONIC SOLID ANGLE 

Note that (see Equations (2-64) , (2-66) and (A-21)) 



*A!I the integrals in this subsection are defined within the theory of distributions 
or generalized functions. Therefore, the concept of the finite part of the integral 
is considered whenever applicable: with this convention the derivative dH/dN c 
is identically equal to zero. 
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Therefore, the analysis for the value of E* on the surface is facilitated by intro- 
ducing the concept of supersonic solid angle defined by 


dQ' = 


R 


2H dQ = - 


2H 


dSr 


3N 


Note that 


f£\ 2 : 
,R / 


(X-X.f 


- 1 


(Y-Y*) 2 + (Z-Z*) 2 


■i-l 


(x-x*r 


(Y-Y*) 2 + (Z -Z*) 2 


+ 1 


Therefore, (R/R 1 )^ H depends only upon (X - X*)/''! (Y - Y*)^ + (Z - Z*)2] , 
Hence, dQ 1 is independent of the distance. This implies that the integral 


= = 2HdC! 




2H 


dE c 


e b sn 


c \ R' 


is independent of the actual surface Eg but depends only upon the topological 
relationship between P* and Eg. In particular, it is apparent that 


0 * = 0 


P* outside Eg 


If P* is a regular point (as defined in Subsection A. 1.1) of the surface 
Eg, the value for Q* may be obtained by replacing the surface Eg with the surface 
Eg, as indicated in Figure A-2. For notational simplicity the point P* is assumed 
to coincide with the origin and the normal N to Eg is assumed to be normal to the 


(A -32) 


(A -33) 


(A -34) 


(A-35) 
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Z-axis. (If this is not the case, the analysis can be repeated with a local coordinate 


system satisfying the above requirements.) With these assumptions one obtains (see 


Figure A-3) 


--0 


3 / 2H 


Z B SN 


c R 


2 2 
Xj - Z A 


= - 2 


5 / 1 


-l x ii j /3“X ax i l R i 


-7 x i- r 


where Rj = (X 2 - Y 2 - Z 2 ) 1//2 . 

3/2 

Consider the rule for the integral of the distribution u / 


f — 

J 0 3/2 
u ' 


du = I 


Ik. f J£_ . 2fW 

e -*0 \ u 3 / 2 e'/ 2 


(A -36) 


(A— 37) 


This yields, for |g'(0)|< co , 


J° -L(_sJ du= j’L -1 g) u- 3 / 2 

Bu V u l/2/ J 0 \ 2 / 


= lim | g 


u- 1 / 2 - (2 ug' - g) 


e -*0 
o 


e~ 1//2 = g(a)//cT (A -38) 


u =0 


It is convenient to rewrite Equation (A-36) as 


Q* = 4 Urn I' 


(A-39) 


1- / 


r 2 9 

x r z x, 


-/ x 2 -, 2 
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or, using Equation (A -38) to perform the integration with respect to Y, 




(A -41) 


Combining Equations (A -39) and (A -41) one obtains 


Q* = 2tt P* on Eg (regular point) 


(A -42) 


(It may be worth noting that if e is replaced with zero, the integrand Is equal to zero 
except at Z = |X||. This indicates that as e goes to zero the integrand in Equation 
(A-41) tends to tt/ 2 times a Dirac delta function located at Z = - JX^ J ») 

It is important to note that the value of Q*, besides being independent of 
Xj (which is to be expected), is also independent of Zj, i.e., the "angle of attack", 
a - tan~* (Z^/X^), of the tangent plane, as long as Z^ < |Xj| . 

The case Z^ > |Xj corresponds to the case of P* inside Eg, for in this 
case the domain of integration in Equation (A-36) is given by the whole area inside the 
circle in Figure A-3. In other words, if P* is inside Eg, Q* is given by Equation 
(A-36) with Zj replaced by } X^ | . Accordingly f is given by Equation (A-40) with 
X^ - e . This yields I e = tt and 
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Figure A-3. Domain of Integration for Supersonic 
Solid Angle for p* on Eg . - 
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(A -43) 


Q* = 4rr P* inside Eg 

(Note the presence of a second Dirac delta function at Z = | ) .) 

In summary 




= 0 

P* outside Eg 



= 2rr 

P* on Eg (regular point) 



= 4tt 

P* inside Eg 

(A-44) 

A.3.2 

VALUE OF E* FOR SUPERSONIC 



Note that 




4rr E* + Q* 

= 4rr 

(A-4 5) 

for P* 

inside and outside E R ; 

therefore. Equation (A-29) may be rewritten as 



4tt$* 




(A-46) 


Repeating the considerations presented in Subsection A. 1.3 one obtains 
that Equation (A -29) is valid for P* outside, inside and on Eg (for both regular and 
non -regular points) if E* is defined as 
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e* = i - 


(A-47) 


In particular, 

E* = — P* on E r (regular point) (A-48) 

2 ' 

is the value used in SOUSSA-P. 

The extension to unsteady flow is easily obtained using the same procedure 
used in Subsection A. 2 for subsonic flow. 
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APPENDIX B 


LAPLACE TRANSFORM IN UNSTEADY AERODYNAMICS 

In this Appendix, issues related to the use of Laplace transform are examined. 
In particular, the question of the truncation of the wake (its rationale and its advan- 
tages) are examined in detail. Finally, the relationship between the Laplace-transform 
analysis for unsteady flow and the complex exponential -analysis for oscillatory flow is 
examined; it is shown that Laplace transform results for Real (p) > 0 may be used to 
study divergent oscillatory flow, whereas convergent oscillatory flows (of the type e^* 
with - oo < t < oo ) are "physically impossible" (the Laplace transform for Real (p) < 0 
exists only for truncated-wake analysis). 

B.l LAPLACE TRANSFORM 

Consider a function f(t) such that 

f(t) =0 (t < 0) (B-l) 

and otherwise arbitrary. The Laplace transform of f(t) is defined by 

f(s) = Lffl = f°° fWe^dt (s € Dp) (B— 2) 

J 0 

where s is the Laplace parameter (complex frequency) and is the domain of 
convergence of the integral. 

Note that f($) is a function of a complex variable: the definition of this 
function outside the domain D^ is obtained by analytic continuation. However, the 
Laplace transform of f(t) is only defined in the domain D^. 
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The inverse Laplace transform yields 


, ~ 1 rot, + i ao ~ , 

f(t) = L~'(f) = 1 f(s)e st ds (B— 3) 

2 * J • 

tt i cr ^ - i co 

where the path of integration is in the domain . In order to evaluate this integral, 
it is convenient to use Jordan's lemma (Reference 24, p. 81): given a single valued 
function such that |g(Re )|-»0, as R-»oo, then 


* ,m f g(s) e st ds = 0 (B-4) 

R -> oo r 

where t is a real constant and T is a semi-circular path with s = Re , where 
-tt/2 < 0 < tt/2 for t< 0 and tt/2 < 0 < 3tt/2 for t> 0. Thus, the path of integra- 
tion may be closed on the left-hand side of the plane for t> 0 (right-hand side for t < 0). 

r+J 

It may be worth noting that if f(s) does not have singularities for Equations 

(B-3) and (B-4) yield immediately f(t) = 0 for t< 0 in agreement with Equation (B-l). 


For nondimensional time T = U t / & , one obtains 

CD 


F(p) = f" HV e-F 1 dT 


(B-5) 


where 


P 


s je/U 


00 


(B-6) 


is the complex reduced frequency (nondimensional Laplace parameter). 

B.2 TRUNCATION OF WAKE 

As mentioned in Section 3, the wake is truncated at a finite distance from 
the body. This point is examined in detail in this Subsection. 
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B.2.1 LAPLACE TRANSFORM OF TRUNCATED FUNCTIONS 


Consider the "truncated functions" (see Figure B-1) 


and 


f,(t) = f(t) 
F,(t) = 0 


yo = o 

f 2 w = f(t) 


t < t 


t > t 


t < t 


t > t 


(B-7) 


(B— 8) 


where t Q > 0. Note that f-j(t) = f 2 (t) = 0 for t < 0. The transforms of f^ and f 2 
are given by 


f,(s) = f'° fWe^'dt 

I Jq 


~ /* 00 

f 2 (s) = J f(0 


-st 

e dt 


(B-9) 


Note that 


f(s) = f,(s) + f 2 (s) 


(B— 10) 


However, if one is interested only in the interval from 0 to t then replacing f with 
f.j yields no difference in the response (since f 2 contributes only for t > t Q ). 

In order to clarify this point further, consider, as an illustration, the 

function 

t > 0 

t < 0 (B-1 1) 


f(t) = e' 
= 0 
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Figure B-l . Truncated Functions. 


The Laplace transform is 


f = f e e" 5 * dt = [ Real (s) > a ] 

*a) s + a 


(B— 12) 


Next consider 


= no 


t < t 


= 0 t > t 


(B-13) 


The Laplace transform of fj is* 


W„ 


t . . 1 / -(s+a)t 

o -crt -st lt _ |, 

e e dt = 1 1 - 

s + cr 


= f - f„ 


(B-14) 


where 


f = 


f 2 ~ 


1 


s + a 

-1 “( s+£y ) t c 


s + a 


[ Real (s) > a ] 
I Real (s) > a ] 


Using Equations (B— 3) and (B-4) and the residue theorem yields 
1 /.a, +ioo 1 


l ” 1 (0 = — / 1 

-t 


e s * ds = 0 t < 0 


and, using t-t instead of t in Equation (B-4), 


L -,y . — l r: 


-(s+a)t 


(B— 15) 


2 rr i * , ar 1 -ioo s+a t 

1 = e at t > 0 (B-16) 


2tti o’^-ico s + a? 


e s * ds = 0 t < t 


= e" at t>f (B-17) 


*Note that the integral in Equation (B-14) exists for all finite values of s. Hence, D p 
coincides with the finite s-plane (f^ is not singular at s = -a). 
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-oo < t < 0 


i.e., 


L“‘(f 1 ) = 0 


= e 


at 


Octet 


= 0 


t < t < co 
o 


(B— 1 8) 


It is apparent, comparing Equations (B-16) and (B— 1 8), that the use of fj 
instead of f yields wrong results only for t > t Q . Therefore, if f^ is simpler to use 
than f, it is legitimate (and convenient) to use fj instead of f in order to obtain the 
response for t < t . 


B.2.2 RATIONALE FOR WAKE TRUNCATION 

In order to understand better the rationale for the truncation of the wake 
it is convenient to consider the Laplace transform of the subsonic integral equation (as 
mentioned in Subsection 3.2, the problem of the truncation of the wake is limited to 
the subsonic flow). For the sake of clarity a simple case is considered here. Consider 
a wing (with straight trailing edge located on the Y-axis) in incompressible flow. 

Then Equation (2-32) for M = 0, yields 0 = 0=0 and Equation (2-52) reduces to 


where 


4rrE(P*)f (P*, T) 



R 



ds B 



dY 



A§ 



dX 


(B— 19) 


(B-20) 
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with 


R = R 
o 


Z = 0 


Using Equations (2-53) and (2-54) with p 
may be rewritten as 


(B— 21) 

1, i.e., II = X. Equation (B-20) 


J W ^D = f A*(Y,T) — {-) 

w *N U \rJ 


dX 


(B— 22) 


The upper limit on the integral is T instead of oo because, according to Equation (2-45) 
A§(Y, T) = A§ te (Y, T-X) (B-23) 

is equal to zero for X > T . The Laplace transform of Equation (B-19) is 


4ttE*5. 


*(P) = “ $ 


B L 

r b/2 ~ 

+ J w dY 

J -b/2 W 

where, using the convolution-integral formula 


- 

~ 1 

~ 9 


Y — 

- $ 

R 

dN 

i R / j 


ds 


B 


L (fj F(T - X) G(X) dX j = L(F)L(G) 


one obtains, from Equations (B-22) and (B-23) 


: W (Y) = Ai; TE r 


3 n 


'0 5N„ \R oJ 


e'P X dX 


(B-24) 


(B-25) 


(B-26) 
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Note that if the wake is truncated at X = T q , Equation (B-20) yields 



which yields, using Equation (B-23) 


(B— 27) 


J$ = J T ° A? te e’P X _JL l±\ dX (B— 28) 

VV J 0 BN,, \ R / 

It is apparent that if A § j E (X T-X) =0 for X > T, is equal to 

for T < T q . Therefore, according to the analysis presented in Subsection B.2.1, the 

truncated wake integral equation may be used for T < T q without any effect on the 

results. This is the rationale for introducing the truncation of the wake; the Laplace 

transform of the solution will depend upon the value of T , but the solution in the 

time domain is identical to the infinite-wake solution for T< T . 

o 

It may be worth noting that J^-» J^y as T o ”*co indicating that the infinite- 
wake analysis may be obtained, in the limit, from the truncated-wake analysis. In this 
case however, the wake integral is expressed in terms of a multivalued function (see 
next subsection). 

B.2.3 ADVANTAGES OF WAKE TRUNCATION 


The integral over the wake is evaluated in Reference 3 for the case 
p = ik. The procedure is briefly outlined here for the more general case p^ = 

Real (p) / 0, in order to show that the integral exists only for p^ > 0 and also 
that the function obtained from the analytic continuation of the Laplace transform 
is a multivalued function which requires the use of branch cuts for the evaluation of 
the inverse Laplace transform. It is also indicated that these problems are eliminated 
by the use of the truncated wake which therefore has the following advantages: 


B -8 



first, the Laplace transform of the integral equation exists in the whole plane (therefore, 
analytic continuation is not required); and second, it is a single valued function 

Consider the wake integral in Equation (2-52) 




d j. 

SN 0 R> ^ 


(B-29) 


For simplicity, assume, in line with the small perturbation assumption, that 
the trajectory of the points of the wake are composed of straight lines parallel to the 
X axis. Then, using Equation (2-53), Equation (B-29) yields 


where 


W 


= / dY f 
TE X 


TE 


a$ te (y, t - e - n) 


3N y \ R 


- Ai TE (Y, T - © - n) 


3 ® 1 


3Ng R 


dX 


The Laplace transform of Equation (B— 30) is 

T w - L ax te t ; 


TE 


■w = J 

A# 


f 00 -p( 0 +n) 

a j 


1 a® l 

J 

X TE 

3Ng I 


! p 

' aNg R 


dX 


(B-30) 


(B— 31) 


(B-32) 


or using for Ng the unit vector in the direction of the positive Z axis and Equation (2-32) 


w 


r°° - P (®+n) 

9 /]) 

M 

- p — 

aR i 

1 e 
X TE 

_ dZ \R) 

P 

9Z R 


dX 


(B— 33) 


B -9 



Note that according to Equations (2-31) and (2-54) 


0 + n 


2 

p 



-X*) + MR + p 2 




(X -X* + MR) + p (X*-X TE ) 


and that 


d e -p(0+n)_ _ e “p(® + n) 
dZ p b z 


(B-34) 


(B— 35) 


Hence, Equation (B-33) may be rewritten as 





“PP(X* g 

dZ 



-£ (X-X* + MR) 
e P — 

R 


dX 


(B-36) 


This integral is of the same type as the one given in Equation (D-l 1) of Reference 3. 


Using the classical transformation 

1 1 _ A + MR 
K o 


(B-37) 


where A = X - X* and 



(Y-YJ 2 + (Z-Z*) 2 



(B-38) 


and noting that 


vw=_l r, 

f"V, L 


p 2 R 2 + A 2 + 2AMR + M 2 (A 2 + R 2 ) 


r 


(B-39) 
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or 


and 


/ 


1 + IT = 


PR 


R 2 + A 2 + 2 AMR + M 2 A 2 

o 


l '/ 2 


_L (MA + R) 
P R o 


-_l (i + M t\ -i yrr? 

dX pR 0 \ R / R 


one obtains 


<B-40) 


(B— 41) 


where 



(X-X. + MR) 

e ° 



l 


-pR U 
oo r o 
e 


dU 


U TE 7l + U 2 


(B-42) 



X -X* + MR 


PR 


X = x 


TE 


(B-43) 


Next, consider the contour indicated in Figure B-2, and note that 




U 


2 


dU 


= 0 


(B-44) 


since the integrand is analytic inside C. Note that using Jordan's lemma. Equation (B-^), 
if the radius of C 2 goes to infinity, the contribution of goes to zero. Also, if the 
radius of goes to zero, the contribution of goes to zero. Also, note that in 
the limit, the contribution of Cg tends to (using V = iU = cosh cr) 
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-ikR U 

_e 

J\ + u 2 


dU 



•kR V 
o 



dV 



-k R cosh a 

e ° do = -K (kR ) 

o' o' 


(B-45) 


where K is the modified Bessel function of second kind of order zero (see Reference 
o ' 

25). Therefore, Equation (B-44) may be rewritten as 




-ikR U 
e °_ 

yr+ u 2 


dU + K (kR ) 
o v o' 


(B-46) 


Combining Equations (B— 36), (B-42) and (B-46) and noting that (Reference 25) 


(where Kj 
for p = ik. 


d K (a) 

— = -«,(*) (B-47) 

da? 


is the modified Bessel function of second kind of first order) one obtains, 


r. -ikp ( x *" x te) 
w = e 


_a_ 

sz 


r - i " ikR o U 

J - e , 

U TE J] + U 2 


dU 


-K^R^ k 


bZ 


(B-48) 


The first term in the bracket can be easily evaluated numerically or analytically (an 
analytical expression is obtained in Reference 3). However, (Reference 25) 
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Kj(o) = |y + + 


oo 

T £ 


2 rn^O m! (m + l)! 


[1 + .1 + 

2 


1 1 


2 m + 1 


'Ck'l 

2 


\2m 


. + 


m + 1 


(B— 49) 


(where y is the Euler constant) is a multivalued function. 


Furthermore, using the same contour of integration, one obtains (using V = i U) 


only for 



(B— 50) 


Real p > 0 (B-51) 

If Real p < 0, Jordan's lemma is not applicable. The contribution of C 2 
goes to infinity (while the integrals over C 3 , C 4 and C 5 remain finite) indicating that 
the integral over Cj does not exist. 

Finally, it is worth noting that if the wake is truncated, then 1^ is given 
by (see Equations (B-36) and (B-42)) 
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•w 


= e 


-PP(X„-X TE ) 3 j U o e -P R 0 U 


sz u. 


TE 


J 1+U 2 


dU 


- P P(X*-X TE ) oo ^ / rU o yn \ 

n = 0 SZ \° U TE 7l+U 2 / n! 


(B-52) 


where the integrals may be easily evaluated by the procedure used in Appendix D 
of Reference 3. (The summation sign and the integral sign may be interchanged because 
the exponential series is uniformly convergent in the interval of integration; this would 
not be legitimate for the infinite wake.) 

Equation (B-52) indicates that 1^ exists for any value of p, since the integral 
of a uniformly convergent series is a uniformly convergent series. 

In conclusion I ^ is a single-valued analytic function of p in the whole p 
plane for the truncated wake, whereas it exists only for Real p > 0 for the infinite 
wake. In this last case the use of analytic continuation is necessary to obtain the 
function in the whole p-plane and this function is multi-valued. 


B.3 OSCILLATORY AND STEADY FLOW 

In this section is shown the relationship between Laplace-transform analysis 

I 

and oscillatory-flow (in particular steady-state) analysis. Consider Equation (3-9) and 


assume Y. and $ . to be of the type 


Y.(T) = Y.e 


ik T 


(B-53) 


and 


§j(T) = i e 


ik T 


(B-54) 
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where 

k = ioX/Uoo (B-55) 

In particular if k = 0 one obtains the steady-state analysis. 

Combining with Equation (3-9) one obtains an expression similar to 
Equation (3-12) with Yj, §j and p replaced by Y j , §. and ik. Therefore, the 
formulation presented here for transient analysis ($. = 0 for T < 0) in the Laplace 
domain is valid for oscillatory flow as well (-00 ^ T <oo) , with the amplitudes of 
the oscillatory solution given by 


Y. = Y. 


§ i = f i 


p =ik 


p = ik 


(B -56) 


This is a well-known relationship between the two types of analysis (the 
same results hold for supersonic flow; see Equation (3-15)): as pointed out in References 
26 and 27, any code developed for oscillatory flow (imaginary analysis) can be used for 
Lap lace -domain analysis (with zero initial conditions) simply by replacing k with 
-ip . 


However, in Equations (B— 53) and (B-54) 


-cd < T < oo (B-57) 

Therefore, the truncation of the wake (which is based upon Equation (3-11); see also 
Subsection B.2), is not applicable in this case. Hence, for oscillatory flows Equation 
(3-13) must contain the complete wake. 

More subtle is the relationship for exponentially growing or exponentially 
damped oscillations; in this case, setting 
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(B-58) 


Y. (T) = Y. e pT 

I I 

f.(T) = §. e pT 


and combining with Equation (3-9) yields Equation (3-12). This yields no problem for 
divergent oscillations, i.e,. 


Real (p) > 0 


(B-59) 


“P n n 

since the terms e -* 0 as II n = p (X n - X^) -* oo (note that II n 
distance of the wake element, 2 , from the trailing edge). 

However, for convergent oscillations, i.e.. 


increases with the 


Real (p) < 0 


(B-60) 


one obtains 



-» co 


as II -> oo 
n 


(B-61) 


i.e., the solution does not necessarily exist. (The results of Subsection B.2.3 indicate 
that indeed the Laplace transform of the wake integral does not exist for Real (p) < 0.) 

It should be noted that this phenomenon is not in contradiction with the physi- 
cal model; for, if the amplitude of Yj at t = 0 is finite, then 

Y. (-co) = co (B-62) 


If Yj (-co) has a "finite effect" at t = 0, then the solution is infinite. 
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B.4 


INITIAL CONDITIONS 


Note that for any arbitrary function f(t) 


L(f ) = sf - f (0) 


(B-63) 


In relationship to Equation (B-63), the question of what restrictions have 
been used on the initial conditions is examined here. 

The flow is assumed to be in steady state for time t < 0. This implies, 
(see Equation 2-45)) 


= 0 (t < 0) 


(B— 64) 


and 


y =0 (t < 0) 
At time t =0 it is assumed 


(B -65) 


q n (o) = 0 


(B-66) 


in order to obtain Equation (4-2 7) from Equation (4-26). 

It is not necessary to assume however that q m (0) = 0. This implies, from 
Equation (4-26), that in general 


yo) / o 


(B-67) 


Equations (B-65) and (B-67) imply, from Equation (2-48) that 


yo) = o 


(B— 68) 
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if and only if ® / 0 (except at P = P*), i.e., for M / 0 (see Equation (2-31)). 
However, in general 


yo) / 0 


(B— 69) 


for M = 0. 

Finally, consider the effect of the initial conditions in obtaining Equation 
(3-12) from Equation (3-9). If M / 0, Equation (B-68) applies and therefore Equation 
(3-12) is correct, since the last term in Equation (B-63) is equal to zero. On the 
other hand if M = 0 then 0 = 0 (see Equation 2-32)) and therefore Equation (3-12) 
is still correct, since in this case 

D., = G. = 0 
|h jn 

(see Equations (3-7) and (3-10)). 

It should be noted also that the case in which Equation (B-64) and (B-65) 
are not valid (i.e., when the flow is not in steady state for t < 0) is examined in 
Appendix F of Reference 3, where it is shown that the formulation is still valid if the 
term given by Equation (F-2) of Reference 3 is added to Equation (2-16). 
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APPENDIX C 

INFLUENCE OF MOTION OF SURFACE 


As mentioned above, in deriving the integral equation, the surface of the 
body has been assumed to be fixed ( e = 0 in Equation (1.10) ). The correct derivation 
of the integral equation for time dependent surface is given in Reference 3. Such 
derivation is not needed here with the only exception of a minor point, which is 
important in deriving the boundary conditions. For, consider the normal-wash 
integral in Equation (2-34) 


Y 


= -(ffi[Y] 0 -L ds 

J7 D 




v i • A, x A, 
o 1 2 


© 


1 |A 1 x A 2 1 dH 1 dH 2 (C-l) 


|A] x A 2 | 

where -f* are the convected curvilinear coordinates (moving with the body surface) and 


A = 
a 


SP 




(C-2) 


are the surface base vectors. 


Note that if the surface is time independent then 


x^l® = xA 2 | 


(C-3) 


However, if the surface is time dependent the two factors are different. Since 
Equation (2-34) was obtained under the hypothesis of time-independent surface, it 
is impossible to determine what is the appropriate value of |A^ x A 2 | to be used. 
Therefore , in order to obtain the correct interpretation of Equation (C-l) it is 
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necessary to make use of the results of Ref. 3, where it is shown that the correct 
expression is (Equation (6-40) of Ref. 3) 


*Y 


Im JVh 9 dE e 

R IV®| 


( 04 ) 


where Y^is the norma! wash in the modified expression of i^, whereas £® is the surface 
described by the equation 


S® = S(P, T - 0) =0 


(C-5) 


S 1 , H 2 ) = P(e\ E 2 , T-0) 


(C-6) 


with 0 given by Equation (2-31). 

H 0 0 

Note that the base vectors of the surface £ are given by 

A® = ^ = TaW_ Fv 0 e- a®] 

a ot o a J 


(07) 


T-0 


(where , P = 9P/6T and 7 © • A® = -^2- ) or , solving for A® ,* 

A® = [*« - P 7 0 • A^ (1 + P • v ©f 1 1 (C-8) 

<* L a ° <* ° Jj-® 

This yields (using P as third base vector, A^ ) 

A® x A® = Ja^ x A^ - (1 + P • V o ©f 1 ^A^ x P %/- A^ x P % / j 

= (1 + P . V^)" 1 UW X A^ (1 + P *V o 0) - A^ X P 0/ + X P 0y ) 

L ° v ° 2 'VJt-0 

*Note that, if a = b + cd • a, then d • a = d * b/(l - c • d) and 
a = b + c d • b/(l - c • d) . 


C -2 



= [(1 + P • V9)" 1 ( X + F • A^ x A<\ @) ] 

= |(1 + P -^e)" 1 |A^ X A^l (N (t) + F • N W V Q 0) ] 


(C-9) 


(where 0y = 7 Q 0 • A^). Note that (see Equation (4-14)) 


a 


iy®i_ 

l».S| 9 ' 


3S 1 


V s 

o 


V SI st |v s 
o ' o 


N^ + F • N^V 0 


V Q 0 


T-0 


(C-10) 


T-e 


Hence 


l vi 


© 


o ' j ^.0 _ 


I v I 


-ds 


+ F • N (f) V 0 

O 




-1 

T-0 

2 


A® x A® 


.«1 J -.2 

da da 


1 + P *7 © 
o 


© 


(C-ll) 


Note that in the common case of a zero -thickness wing, at zero angle of 

attack, with motion in the direction of the normal, 7 0 is in the plane of the wing, 

• * 

whereas P is along the normal N and hence P • 7^0 = 0. In general for small -amplitude 
oscillations at reduced frequencies of order one, one obtains 

1 P - 7 o ©| < < 1 (C-12) 

This approximation is used in this report. 

Combining Equations (C-4), (C-ll) and (C-12) and noting that, for small- 
amplitude motion, R®z R, one obtains 
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or 




d H' 



(C-13) 


(C-l 41 


where R, |Aj x A 2 1 and ds are evaluated at time T or, more conveniently in the 
undisturbed configuration. 


The results of this subsection may be summarized by saying that Equation 
(C-l) is still valid if Y 0 is defined as 


Y 


© 


1A 1 x A 2 |® 

|Aj xA 2 



(C-l 5) 


where Y^has been introduced in Equation (C-4).. Note that Equation (2-44) 
applies to Y^ as well, and therefore 


(C-16) 


Note also that 


|A 1 x A 2 | 0 _ |a ] x a 2 | 0 

|A 1 xA 2 j 


(c-l 7) 


Therefore in order to take into account the motion of the surface it is convenient to 
retain Equation (C-l) with Y given by 


C-4 



1 


(C-l 8) 


V - 1 |, 

U oo 


i°S^ x 4° I 


M' 


l a i X a 2 | 


In Equation (C-14) the superscript (t) is used to indicate evaluation with time- 
dependent surface (e r 0) ' n Equation (1-10)). The base vectors a^ are evaluated 
from the time-independent surface (e = 0 in Equation (1-10) ). 
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APPENDIX D 


BODY-AXIS BOUNDARY CONDITIONS 


In order to extend the formulation of Section 4 to the body-axis formulation 
(which is used in flight dynamics, see for instance Reference 28), it is necessary to 
express the downwash in terms of the generalized coordinates and generalized velocities. 
This relationship can be obtained by using the boundary conditions as obtained by 
combining Equations (4-16) and (4-18) to obtain 



(n + An) 


(D-l) 


or 


Y 


Y c + Y 


U 


(0-2) 


where 



(D-3) 


(D-4) 

It should be noted that in the body-axis formulation u is the displacement 
(of a point P on the surface of the body) evaluated with respect to the frame of 
reference ig, jg, kg. Therefore, including the motion of the body axes to Equation (4-13) 
with u given by Equation (4-23) one obtains (see Reference 28) 


and, neglecting higher order terms 

V. 


Y u = 


B 


n - ! • An 


U 


CD 
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V B 


= U 


CD 


6 

E 

n = 1 


v M 
n n 


(R) 


if) 


+ i 


M 

E 


q M (§' 




( 0 - 5 ) 


m 


where typically M are the nondimensional natural modes of vibration and a are the 
corresponding nondimensional generalized coordinates, whereas v n are the nondimen- 
sional components of the translational and angular velocities of the body-axis frame of 
reference whereas /Vv^' are the nondimensional six rigid-body mode shapes. Note that 
the nondimensional translational velocities are obtained from the translational velocities 
by dividing by (the corresponding modes are already nondimensional). On the 
other hand, the nondimensional angular velocities are obtained from the angular veloci- 
ties by dividing by U /i (the corresponding modes are made nondimensional by dividing 
by i). 

Combining Equations (D-4), (D-5), and (4-24) yields 


T u - E 

n = 1 


v M 
n n 


(R) - 

v ' • n 


M 
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m 
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— • 

n H 


- in 
m m 


1 q, 
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CD 


The Laplace transform of Equation (D-6) is* 


¥ = 


6 

E 

n = 1 


M< R > • n 
n 


M 

E 

n = 1 


[p M, 


n - in 


m 


(d-7) 


Equation (D-7) may be rewritten as (the subscript h indicates the evaluation at the 
center of the element £^) 

tv - [ E £n ,V) ] (V + [ E ? m ’ 0) ] (V (D- 8 ) 

*|t is assumed q m (0) ~ 0; see Section B.4. 
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where 



and 


£0,Q) 

hm 


[p M * n - An 
m m 



(D-9) 


(D— 10) 
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APPENDIX E 


CLOSED-WAKE PHENOMENON 

Consider, for instance, the wake emanating from the perimeter of the base 
of a truncated cylinder. This wake separates the flow field into two parts; one inside 
the wake and one outside the wake. This type of wake will be referred to as a "closed 
wake". In this Appendix a few mathematical details relative to the analysis of a 
closed-wake configuration are discussed. It will be shown that in this case the integral 
equation (Equation (2-52) is singular (i.e., the homogeneous integral equation has a 
nontrivial solution, and hence the nonhomogeneous integral equation has a nonunique 
solution). For simplicity only the steady incompressible flow is considered here. 

The integral representation of the potential for steady incompressible flow is 
given by (see Equation (2-52)) 



If p* is on Eg, the function E(p) assumes the value E* = 1/2 and Equation 
(E-l) reduces to an integral equation for the potential on the surface of the body. 

As previously mentioned, whenever the wake separates the flow field info 
two separate regions, one inside and one outside the wake, the integral equation is 
singular. This problem occurs in the formulation for bodies with a blunt back such as a 
projectile or a fuselage with truncated base. A similar problem occurs in the formulation 
for building aerodynamics and actuator -disk aerodynamics. 
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In order to show that the integral equation for closed-wake configurations 
is singular, it is convenient to obtain Equation (E-1) a different way. Consider Figure 
Figure E-l(c), and note that the value E = 1 outside the body may be obtained as 
the sum of the function E. for inner flow (i.e., inside the wake) and the function E q 
for the outer flow (i.e., outside the wake and the body) or 


E(p) = E.(p) + E o (p) 


(E-2) 


Applying Green's theorem to the regions shown in Figures E-l(a) and E-l(b) 


yields 


and 


4tt E. (p*) cp(p*) = 


1 


(1 . 9 A 

dn. ^r J dn. I r 


dE. 




(E— 3) 


(E-4) 


Note that 


n. = n ; n = n on £ 


B 


n. = -n ; n = n on E 
i u' o u 


W 


(E-5) 


Thus, adding Equations (E-3) and (E-4) and using Equations (E-2) and (E-5) yields 


4nE(p*)cp(p*) = "fi— — (-) 

[_dn \r j bn \r f m 


dE 


// K~Vi> 


± i\ dE 
9n u\ r 


(E-6) 
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Figure E-l. Geometry for Green's Theorem for Closed-Wake Configuration. 
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Note that Equation (E-6) coincides with Equation (E — 1) with Acp = cp Q - cp. . In order 
to show that the operator in Equation (E-6) is singular, it is sufficient to show that the 
homogeneous equation (i.e., with dcp/dn =0) 


4tt 



has a nontrivial solution. , In order to show this, consider the function 


(E~7) 


cp (p) = 0 p outside £. 

NT ' 

= 1 p inside E. (E-8) 


If Scp/dn = 0 (homogeneous problem), this function is a trival solution 
(i.e., 0 = 0) of Equation (E-4), whereas for Equation (E-3) one obtains (note that 

E i cp NT = E i^ 

4 " E >^^(;) dE > (E - 9) 
y I \ / 


which is identically satisfied, since (indicating with 0 the solid angle) 

\ d£. s n = 0 outside 2. 

I ' 

= 2n on E. 

= 4n inside E. (E-10) 

Since iso nontrivial solution for both Equations (E-3) and (E-4) with dcp/dn = 0, 
it is a nontrivial solution for their sum. Equation (E— 7). Hence, the operator in 
Equation (E-1) is singular; therefore, if cp satisfies Equation (E-l), then 
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cp = 9 + Cep 


(E-ll) 


NT 

(where C is an arbitrary constant) also satisfies Equation (E-l). 

Finally note that in the numerical formulation. Equation (E-l) is replaced by 
Equation (3-12) with p = 0,. However, the doublet integral represents solid angles and the 
solid angle is evaluated exactly with the use of the hyperboloidal elements. Thus the 
discrete form of Equation (E-10) is still valid exactly. Therefore, even the discrete 
system is singular, i.e., the determinant of the system given by Equation (3-12) with 
p = 0 is equal to zero. This implies that the vector 

{^i } = { 9 NT^V} (£~^) 

is a nontrivial solution for Equation (3-12) with right-hand-side equal to zero, or that 
if 9. is a solution to Equation (3-12), then 

{ 9 i} = { 9 i} + C I 9 !'' 1 } ( E “ 13 ) 

(where C is an arbitrary constant) is also a solution to Equation (3-12) (i.e., if the 
solution exists, it is not unique). 
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